arXiv: 1501.04515v2 [cond-mat.soft] 27 May 2015 


On the relevance of disorder in athermal amorphous materials under shear 


Elisabeth Agoritsas, 1,2, * Eric Bertin, 1 ’ 2 Kirsten Martens, 1 ’ 2 and Jean-Louis Barrat 1,2 

1 Univ. Grenoble Alpes, LIPHY, F-38000 Grenoble, France 
2 CNRS, LIPHY, F-38000 Grenoble, France 
(Dated: May 28, 2015) 

We show that, at least at a mean-field level, the effect of structural disorder in sheared amorphous 
media is very dissimilar depending on the thermal or athermal nature of their underlying dynamics. 
We first introduce a toy model, including explicitly two types of noise (thermal versus athermal). 
Within this interpretation framework, we argue that mean-field athermal dynamics can be accounted 
for by the so-called Hebraud-Lequeux (HL) model, in which the mechanical noise stems explicitly 
from the plastic activity in the sheared medium. Then, we show that the inclusion of structural 
disorder, by means of a distribution of yield energy barriers, has no qualitative effect in the HL 
model, while such a disorder is known to be one of the key ingredients leading kinematically to a 
finite macroscopic yield stress in other mean-field descriptions, such as the Soft-Glassy-Rheology 
model. We conclude that the statistical mechanisms at play in the emergence of a macroscopic yield 
stress, and a complex stationary dynamics at low shear rate, are different in thermal and athermal 
amorphous systems. 


I. INTRODUCTION 

Understanding the nature of the plastic response of 
amorphous media to externally applied forces, and their 
resulting mechanical and rheological properties, is a chal¬ 
lenging issue in current material science research. Inter¬ 
estingly, materials that seem at first sight very differ¬ 
ent, like amorphous solids ( e.g . metallic glasses, polymer 
glasses, granular materials) or yield stress fluids {e.g. gels, 
foams, dense emulsions), share qualitatively a very sim¬ 
ilar yielding behavior in their plastic response to shear. 
Such a resemblance between the dynamics of hard and 
soft materials has already been highlighted in the field 
of crystal studies, where L. Bragg and J. F. Nye used 
monodisperse bubble rafts in compression experiments, 
in order to mimic the dynamical properties of crystalline 
structures [1]. 

For dense, structurally disordered materials, it has 
been known for a long time that the yielding process 
starts at a well-defined material-dependent stress thresh¬ 
old. The ensuing flowing regimes are traditionally char¬ 
acterized by constitutive laws that describe plasticity as 
a homogeneous steady flow. It was A. Argon who first 
introduced the idea of localized shear events on a micro¬ 
scopic level as the physical mechanism underlying plastic¬ 
ity [2, 3], similar to the existence of defects in crystalline 
materials. This idea of localized events is also at the 
basis of the Princen theory of foams, that features the 
so-called “Tl” events, small rearrangements of bubbles 
in a disordered foam, that collectively generate the flow 
[4]. These pioneering ideas inspired many works both 
in numerics and in experiments. Taking advantage of 
the modern developments in particle tracking techniques 
and of the boost in computer power, the idea of localized 
plastic events has been widely validated as a microscopic 
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scenario. Important experimental verifications have been 
given e.g. in the field of colloids [5], gels [6] and granular 
materials [7], as well as in simulations of the response to 
shear in glasses at the particle scale [8, 9]. 

The local yielding picture has also led to the develop¬ 
ment of a variety of models at a mesoscopic scale, with 
currently an increased research activity in order to con¬ 
nect qualitatively and quantitatively their predictions, 
both to simulations starting from a microscopic model¬ 
ing, and to the observed macroscopic mechanical and rhe¬ 
ological properties. For a recent review on this topic see 
ref. [10]. The most important challenge for these meso¬ 
scopic approaches is to find the correct way to model the 
local yielding dynamics and to implement the effect of the 
long-range elastic response of the surrounding medium 
to the locally plastic zones. Numerical studies on meso¬ 
scopic models with a spatial resolution do actually im¬ 
plement these interactions explicitly, using the Eshelby 
theory of elastic response to a local deformation [11-13]. 
The resulting mechanical noise is in this way triggered by 
the yielding dynamics itself and is thus, by construction, 
self-consistent. When constructing mean-field descrip¬ 
tions, it is then crucial to develop a mesoscopic picture 
integrating accurately into the evolution equations this 
mechanical noise and the associated activation processes. 
This issue is even more important for the so-called ather¬ 
mal systems, where the stress relaxation due to thermal 
noise can be neglected, and new plastic events are trig¬ 
gered solely by the macroscopic shear and the resulting 
mechanical noise. 

At a phenomenological level, there are several ways to 
describe the mechanical noise. One of the first propos¬ 
als was put forward in the Soft-Glassy-Rheology (SGR) 
model [14, 15], by assuming that the mechanical noise 
acts as an effective activation temperature x, controlling 
an Arrhenius-like rate of plastic events. Combined with 
structural disorder, this description yields a broad dis¬ 
tribution of relaxation time scales for a sufficiently low 
effective temperature. It leads in particular to a complex 
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fluid behavior, to the emergence at sufficiently small x 
of a Herschel-Bulkley type of rheological curve [16]. In 
other words, its macroscopic stress gm as a function of 
the applied shear rate 7 follows a power law with a thresh¬ 
old: g m ( 7 ) (jy-GR yjSGR ^(i-x) ^ p- g exponent de¬ 
pending explicitly on the effective temperature. This ap¬ 
proach has recently been brought into question, at least 
in the scope of strictly athermal dynamics [17], by point¬ 
ing out that the subtle interplay between the different 
time scales at play in the sheared dynamics can actually 
jeopardize the exponential, Arrhenius-like rate of plastic 
events, at the core of the SGR model. 

Another way to describe the mechanical noise, as pro¬ 
posed in the Hebraud-Lequeux (HL) model [18], is to 
model it more explicitly as a diffusion of the local stress 
with a diffusion coefficient proportional to the plastic ac¬ 
tivity, z.e., the rate at which plastic events occur. As we 
shall see below, this modeling of the mechanical noise dif¬ 
fers from the SGR approach not only because the noise 
amplitude is a dynamical variable, but also because the 
noise does not act on the same physical observables in 
the SGR and HL models. This second model predicts in 
particular three different scaling regimes for the rheologi¬ 
cal law gm{ 7 ) controlled by a single coupling parameter, 
including a Herschel-Bulkley behavior of fixed exponent 
1/2. It is important to note another crucial difference be¬ 
tween these models: the SGR model includes the struc¬ 
tural disorder as a key ingredient, without which no com¬ 
plex rheological behavior can emerge, while the original 
version of the HL model does not include any structural 
disorder. 

In this paper, we study a generalization of the HL 
model including a distribution of yield energy barriers, 
and we show that its predictions for the rheological law 
at low shear rate are qualitatively robust with respect 
to the introduction of structural disorder. The paper 
is organized as follows. We first present in sect. II a 
toy model that we use as an interpretation framework in 
order to distinguish between the thermal and athermal 
types of noise in the description of sheared amorphous 
media. Then we recall in sect. Ill the definition and 
main properties of the original HL model, before study¬ 
ing in sect. IV its generalized version in the presence of 
structural disorder. We conclude that the physical mech¬ 
anisms that lead to a non-linear macroscopic response 
to shear are different in thermal and athermal sheared 
amorphous systems. Our mean-field mesoscopic analysis 
suggests that the impact of structural disorder strongly 
depends on the nature of noise (mechanical versus ther¬ 
mal) in the local stress dynamics. 

II. MECHANICAL NOISE VERSUS THERMAL 
DYNAMICS 

In order to understand the distinction between the 
thermal or athermal nature of noise, it is useful to con¬ 
sider the following simple description, distinguishing two 


different local degrees of freedom on which a noise might 
act. In a coarse-grained description, we decompose the 
system into small boxes, for which we can define the local 
strain £ (that we assume to be scalar for simplicity) and 
the local stress g at each time £, under an external shear 
rate y(t). We emphasize the different interpretations of 
these two variables: £ characterizes the local configura¬ 
tion of the system within the box, while g characterizes 
the forces exerted by neighboring boxes, in line with the 
continuum mechanics interpretation of stress. 

If the local element behaves as a linear Hookean spring, 
g and £ are linearly related. However, as plasticity results 
from irreversible rearrangements of particles at a micro¬ 
scopic scale, a complete mesoscopic description should 
also include, in addition to the evolution of the stress, 
the dynamics of an the ’internal’ degree of freedom £. In 
what follows, we first define a toy model for the coupled 
dynamics of £ and cr, and we use it as an interpretation 
framework in order to infer and discuss the assumptions 
made on the underlying dynamics of the strain £ for a 
given effective dynamics of the stress g. 


In the absence of stress, we first assume that the evolu¬ 
tion of the coarse-grained internal variable £ is described 
as a dynamics in a potential Vo {£), characterizing the lo¬ 
cal potential energy landscape of the possible configura¬ 
tions within the box. For simplicity, we also assume this 
dynamics to be overdamped. Secondly, we further as¬ 
sume that in the presence of a local stress cr, the potential 
Vo(£) is modified by a linear contribution proportional to 
< 7 , becoming thus V a (£) = Vq{£) — g£. In other words the 
stress tilts the potential energy landscape. Thirdly, we 
assume that the local strain is subjected to a Gaussian 
white noise £t(£) ~ interpreted as a thermal noise as it 
acts on the local configurational degree of freedom £ - of 
zero mean and fixed variance, satisfying 

(£r(£) £r(0> = 2T(5(t - t') (1) 

with T the ‘temperature’ of the system, ( the damping 
coefficient controlling its relaxation and (...) denoting 
the statistical average over the noise. This term describes 
the thermal contribution to the forces exerted by the sur¬ 
roundings, due e.g. to acoustic-like vibrations. It can be 
significant for a small system at non-zero temperature, 
even in the absence of a macroscopic driving, but would 
be strictly zero in an athermal system. The evolution 
equation for £ in this tilted potential then reads: 

Cd t £ = -+ <r(t) + £r(£) (2) 

In this picture, a fluctuation of £ within a local minimum 
of V a (£) corresponds to an elastic local displacement of 
particles, whereas a jump above an energy barrier to¬ 
wards a new local minimum corresponds to a plastic lo¬ 
cal rearrangement, which can be triggered either by the 
‘thermal’ noise or by the fluctuating local stress cr(t). 
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Figure 1. Summary of the mean-field picture presented in sect. II, centered on two coarse-grained degrees of freedom: the local 
deformation £ and the local stress a (its shear component), whose coupled dynamics is governed by eqs. (2)-(4). Left: Starting 
from a microscopic description of an amorphous material, subjected to an external shear rate 7 (£), we distinguish for a given 
mesoscopic box the local configuration £ of the particles inside the box, from the local stress a resulting from the particles 
outside the box. Center: In this simplified picture, if these two degrees of freedom were considered separately, £ would relax 
in the configuration potential Vo(£) accounting for the local structural disorder, whereas a would fluctuate around a sawtooth 
behavior controlled by a given rate u(cr,£,£) of plastic events. Right: Their dynamics are actually coupled, tilting the bare 
potential Vo with the fluctuating stress, and reciprocally relating the local yield stress a c to the closest inflection point of Vb; 
the thermal noise and the mechanical noise £ p i are defined as the noise of the individual Langevin dynamics of the two 
degrees of freedom £ and a. 


The dynamics of the local stress a(t) is assumed to be 
driven by the externally imposed shear rate which in¬ 
crease the stress, and by the yielding events which relax 
the stress. The dynamics is also affected by distant plas¬ 
tic rearrangements, that we model through a white noise 
£ p i (t) acting on a. We thus postulate the following dy¬ 
namics, 

d t <j = G07W + F(cr, £, £) + £ p i(t) ( 3 ) 

where Go is the shear elastic modulus, £ = d±£, and 
F(cr, £, £) is a coupling term between stress and deforma¬ 
tion, leading to stress relaxation during a yielding event. 
The form of this coupling term is not known, although 
one might guess that it may crucially depend on the de¬ 
formation rate £, which becomes significantly larger after 
£ has crossed the energy barrier. For simplicity, £ p i (t) is 
assumed to be a Gaussian white noise, with a dynamical 
amplitude D p \(t) that may (slowly) depend on time: 

(£ P i (t) £ P i(*')> = d pi( t) ~ t') ( 4 ) 

In the following, we further simplify the dynamics, and 
replace both the explicit dynamics of £ given in eq. (2) 
by an effective stochastic rate v(a) with which the local 
stress (j instantaneously relaxes to zero. We thus end up 
with an effective hybrid stochastic evolution equation for 
< 7 , of the form 

d t cr = Go7(£) +£ p i (t) (5) 

a i-» 0 with rate z/(cr, £, £) (6) 

The first term on the right hand side of eq. (5) describes 


the continuous elastic load of the local stress due to the 
shear rate 7 (t) (Go is the elastic shear modulus), while 
the second term £ p i(t) is a noise accounting for the ef¬ 
fect of distant plastic rearrangements, and as such can 
be interpreted as a mechanical noise. Combining eqs. ( 2 ) 
and ( 6 ), it is already clear (see also discussion below) 
that this noise has a cumulative effect on the dynamics 
of £, very different from the uncorrelated thermal noise 
£t- I n order to fix the diffusion coefficient D p \(t), a clo¬ 
sure relation has to be found, involving for instance the 
rate of plastic events in the system. For example, in the 
HL model [18], a simple closure relation is provided by 
the physically reasonable assumption that the diffusion 
coefficient is proportional to the average rate of plastic 
events. This whole physical picture is summarized in 
fig. 1 . 

If we make a change of variables and use, instead of 
the strain the local energy barrier E (i.c. the distance 
in energy to the yield point, which depends on a and £), 
the dynamics given by the eqs. (5)-( 6 )-(4) leads, even¬ 
tually, to the following evolution equation for the joint 
distribution of the local stress a and local yield energy 
barrier E : 

d t V(E, a,t)=- Goi(t) d a V + D pl (t) d 2 a V 

— v(a, E)V + r(t) 5(a) p(E) ' ^ 

with the plastic activity T(t) defined as the averaged plas¬ 
tic rate (v(<j,E)}^ = J dE f dcr v(a, E)V(E, cr, t), and 
p(E ) the probability distribution of energy barriers. Al¬ 
ternatively, E could be replaced by the local yield stress 
cr c or any other relevant feature of the configurational 
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potential Vo(£). Such an evolution equation, with this 
specific structure, is precisely the starting point of the 
SGR [14, 15] and HL [18] models, two prototypal mean- 
field models for sheared amorphous systems. In our toy 
model picture, the assumptions made on the underlying 
dynamics of £ thus translates into a given set of effective 
parameters {D p \(t), z/(cr, E), p{E)}. 

In the above formulation, the fact that two different 
types of noise may be present appears clearly. On the 
one hand, the noise in eq. (2) may be interpreted 

as a thermal noise, as it act on the local configurational 
degree of freedom L On the other hand, the noise £ p i(£) 
expresses the effect of the mechanical noise, originating 
from distant plastic events. This picture is of course a 
simplified local mean-field approximation, because (i) the 
local deformation £ is only an effective coarse-grained 
variable, which does not describe all the positions of the 
particles in the box considered, (ii) the relative prox¬ 
imity of distant plastic events is encoded in non-trivial 
correlations of the noise £ p i, which are completely ne¬ 
glected via the Gaussian white noise simplification, and 
(Hi) the local stress does not necessarily fully relax to 
zero after a plastic rearrangement, nor does it relax in¬ 
stantaneously. Nevertheless, the important point is that, 
in our picture, both noises act on different degrees of free¬ 
dom, and thus do not play an equivalent role. The noise 
£ T (t) acts on the configurational variable £, which experi¬ 
ences a return force making hard to overcome the energy 
barrier (for a fixed energy landscape), leading to long and 
broadly-distributed escape times through the Arrhenius 
relation assumed in the SGR model. On the contrary, 
the noise £ p i (t) acts on the local stress a without any 
return force, so that reaching the value a c at which the 
potential energy landscape changes of minima is compar¬ 
atively easier, and takes a much shorter time. An alter¬ 
native formulation [17] is to consider the fluctuations of a 
as a mechanical noise acting on the deformation £. How¬ 
ever, in this view, the correlation of the noise does not 
decay to zero on a short time, but rather increases due to 
the persistent deformation induced by plastic events (or 
equivalently, due to the absence of recoil force acting on 
a). An important physical consequence of this property 
is related to the influence of disorder, as we shall see in 
Section IV. 

Within the present framework, the situation without 
mechanical noise (£ p i(t) = 0) corresponds to a purely 
‘thermal’ dynamics of the local strain i. In this case, the 
local stress has a sawtooth behavior without fluctuations, 
whose stress drops can be identified as fast relaxations of 
the strain i into a new local minimum of the potential 
energy landscape. If we approximate a given well of the 
unstressed potential Vq{£) by a harmonic potential cen¬ 
tered on £q, we have in presence of a local stress a that 
V a (£) ~ \k{£ — £ o) 2 — cr(£ — £o) and the local strain £ a 
corresponding to its minimum is given by £ a — £q = a/k. 
If we assume furthermore that plasticity is dominated by 
plastic events triggered by the thermal noise £y(t) (and 
not by the local stress), then the dynamics of £ given by 


eq. (2) can be implicitly reduced to an Arrhenius escape 
rate above an energy barrier A E(a). The effective plastic 
rate in the stress evolution equation (6) can then be writ¬ 
ten as v(AE(a)) oc exp [—A E(a)/T]. Since the stress dif¬ 
fusion coefficient has been assumed to be D p \(t) = 0, this 
Arrhenius-like activation description of plastic events is 
qualitatively similar to the SGR model [14, 15], in which 
the structural disorder encoded in p(E) has been shown 
to play a crucial role on the rheological properties. 

In contrast, the situation described by the HL model 
corresponds to a purely relaxational dynamics, without 
thermal noise, so that £ simply relaxes to the local mini¬ 
mum of the potential V a (£). In this picture, a threshold 
cr c is defined as the maximum value of | dV ^ | of the cur¬ 
rent barrier to overcome, and as such it determines the 
minimum value of local stress allowing for a change of 
potential well -a plastic rearrangement within the meso¬ 
scopic region- in the athermal case, z.e., = 0)- The 

rate v is simply taken as zero if the local stress a is be¬ 
low the threshold a c and takes a constant value 1/r if 
the local stress exceeds this threshold. Moreover, in the 
original HL model [18], a c can take only a single typical 
value, whereas we will consider more generally a distri¬ 
bution of threshold values p(a c ). Note that in the HL 
model, the degree of freedom £ is not explicitly described 
either, but including it allows for a better understand¬ 
ing of the origin of the local rearrangements in terms of 
the evolution of the local potential energy landscape. We 
can for instance infer, at least formally, the distribution 
p(cr c ) from the distribution of the disordered potential 
V [Vo CO] if self. To sum up, in the HL model, only the 
mechanical noise £ p i (t) is taken into account and there 
is no ‘thermal’ noise on the configurational variable £, so 
in our interpretation framework the HL description cor¬ 
responds an athermal dynamics of a sheared amorphous 
system, whereas an Arrhenius-activated plasticity corre¬ 
sponds to a thermal dynamics. The main goal of this 
paper is to investigate the role of disorder, encoded in 
the distribution of barrier heights E or yield stress cr c , 
on the rheological properties in the case of a pure ather¬ 
mal dynamics as described by the HL model. 

Let us conclude this section with a word of caution. 
First, in the simplified framework presented above and 
summarized in fig. 1, it appears rather natural to iden¬ 
tify the temperature T with the ‘real’ temperature, z.e., 
to set it to zero in athermal systems. Nevertheless, most 
interpretations of the SGR model introduce the notion of 
an effective temperature associated with the mechanical 
noise, which has the same role, but may be unrelated to 
the physical temperature. Some other models, such as 
the Shear-Transformation-Zone (STZ) model [19], also 
introduce the notion of an effective temperature, which 
is however not used, in general, to compute an activation 
rate but rather as an internal variable that characterizes 
the state of the material. Secondly, in the view presented 
above, the mechanical noise is only described by the fluc¬ 
tuating term in the stress evolution, and is thus totally 
absent from the equation of the strain evolution (except 
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for its coupling to the fluctuating stress a). Whether or 
not a remnant of the mechanical noise should also be con¬ 
sidered at this level, due to the very strong simplification 
that replaces the complex energy landscape of a few tens 
of particles with a single degree of freedom £, remains an 
open question. 


III. THE HEBRAUD-LEQUEUX MODEL 


first equality in ( 10 ) is imposed in general by the nor¬ 
malization daT(a, t) = 1 of the PDF- and is then 
quantified by the proportion of overstressed regions in the 
system. The Dirac distribution S(a) in ( 8 ) corresponds to 
the full relaxation condition stated in ( 6 ), and this con¬ 
dition could be replaced more generally by a distribution 
A(cr) of the stress after relaxation. Finally, and this is 
the key ingredient of the model, the diffusion coefficient 
(4) is assumed to be proportional to the plastic activity, 
hence the following linear closure relation: 


Although the HL model [18] has been previously stud¬ 
ied in great detail in the mathematical literature [20-26] , 
no concise explicit account of the derivation of its rhe¬ 
ological law is available so far in the physics literature. 
Before generalizing the HL model by including a distri¬ 
bution of threshold stresses (see Sect. IV), we find useful 
to provide the main steps of the derivation for the stan¬ 
dard HL model. In this section, we thus recall first the 
definitions of the HL model, second the main steps of 
the derivation at fixed shear rate 7 (and further in the 
limit 7 —)> 0 ) of its stationary solution, third its corre¬ 
sponding predictions for the macroscopic stress <jm{ 7 ), 
and fourth its connection with the Kinetic-Elasto-Plastic 
(KEP) model [27], since all these points will prove use¬ 
ful for the study of its disordered generalization in the 
next section. We detail in particular the scalings of the 
macroscopic stress, with their associated prefactors, in 
the different limits of interest, expressing them in a way 
that will be systematically generalized in sect. IV. 


A. Definition of the HL model 

Introduced closely after the SGR model, the origi¬ 
nal HL model [18] is defined by the following evolution 
equation for the probability distribution function (PDF) 
V(a,t) of the local stress a at time £, under an external 
shear rate 7 (£): 

d t V(a, t)=- G„7(i) + D Hh (t) d 2 a V 

+ T(t)6(a) 

It corresponds to the hybrid stochastic dynamics de¬ 
fined by eq. (5)-(4) for the mean-field local stress, with 
Dpi (t) = Dhl(^), and the following specific rate z/hl and 
plastic activity T(t): 

i / hl(ct, <t c ) = - 6 ( \a\ - a c ) (9) 

T 

Y(t) m(un L (a,a c )) v(att) * 1 [ da'V(a',t) (10) 

T J\(J'\>(J C 

where 0 is the Heaviside function and S the Dirac dis¬ 
tribution. The choice (9) assumes that there is a single 
typical value of the threshold stress a c and that the rate 
of plastic event can be approximated by a fixed value 
1 /t in any overstressed region. As for the plastic activ¬ 
ity, it is defined as the mean rate of plastic events -the 


Duhit) = ocT(t) ( 11 ) 

where a > 0 is for the time being an ad hoc parameter 
of the model. Since the plastic activity depends on the 
PDF, this relation implies that the evolution equation ( 8 ) 
is nonlinear in V(cr,t), the nonlinearity being encoded 
in the diffusion coefficient and leading to the non-trivial 
features of the HL model. 


B. Stationary solution at fixed shear rate 


We focus exclusively on the case of constant shear rate 
7 , which has also been studied in the mathematical liter¬ 
ature [20, 23-25], proving in particular the existence and 
uniqueness of its stationary solution at a constant shear 
rate, in the specific limit 7 -» 0 . 

Assuming that the stationary PDF V st (a) exists, its 
corresponding plastic activity r st and stationary diffu¬ 
sion coefficient Dhl are well-defined, hence the determi¬ 
nation of the stationary solution of the HL model pro¬ 
ceeds generically in the following way. (i) One does not 
take into account the closure relation eq. ( 11 ), and con¬ 
siders Dhl as a fixed diffusion constant D. (ii) The sta¬ 
tionary solution of eq. ( 8 ) can then be determined, and 
in this specific case one has to solve a second-order dif¬ 
ferential equation with constant coefficients on each of 
the intervals (—oo, —<j c ), (—cr c , 0 ), ( 0 , a c ) and (<j c ,+oo), 
connecting them using the continuity of V st (a) and of its 
derivative at a = =bcr c . (Hi) V s t (cr) has actually an overall 
multiplicative constant, that can be identified with the 
(as yet unknown) plastic activity T st defined by ( 10 ); us¬ 
ing the normalization condition daV st (a) = 1 , one 
ends up with an equation of the form: 


r s t r = — 

/<T 0 


Dt 

(■JDt Yto±L \ 

yv-LJT, Dr j 


( 12 ) 


where f ac (x,y) is a known function (cf. Appendix A). 
(iv) At this stage, the closure relation ( 11 ) can at last 
be taken into account, yielding either D = 0 , or a finite 
diffusion coefficient according to the condition: 

(13) 


from which D = Dhl can be determined uniquely as a 
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function of the shear rate 7 and of the coupling parameter 

Note first that this procedure is quite generic, in the 
sense that another choice for the rate z/(cr, a c ) in eq. (9) 
will only modify the functional form of f Gc , and secondly 
we are free to choose a different closure relation than 
(11) starting from eq. (12). The last equation (13) has 
actually a straightforward geometrical interpretation, il¬ 
lustrated in fig. 2. The two arguments of the function 
f<j c , respectively x = V Dr and y = Go 7 t/x 2 , are natural 
choices given the structure of the HL equations (8)-(9), 
that we thus keep on purpose in all our results. 



Figure 2. Schematic plot of the function j Gc y^/Dr, j 

defined in eq. (12), in the presence of a constant external shear 
rate 7 > 0. According to eq. (13), this plot provides a geomet¬ 
rical interpretation of the solution x = V^hl t imposed by 
the closure relation (11), as the intersection between the func¬ 
tion f ac and the horizontal line a. The black dashed parabola 
is the limiting case at 7 = 0 of eq. (14), whose minimum de¬ 
fines the critical value a c ■ The three curves with increasing 7 
correspond to the more general case of eq. (16). The different 
cases for lim^^o Dhl( 7, a ) listed in eqs. (17)-(19) correspond 
to the intersections with f<j c which occur above a c (area 1), 
below di c (area 2) or exactly at d c . 


C. Stationary solution in the absence of shear rate 

In absence of shear rate (7 = 0), the stationary 
PDF is symmetric so it predicts no macroscopic stress: 
o M = da aP st (a) = 0. This case is nevertheless a 
benchmark of the model, as it exhibits two regimes for 
the diffusion coefficient depending on the value of the 
coupling parameter a in the relation (11). Indeed, the 
HL model at 7 = 0 yields: 

2 

x = \/Dt, f ac (x, 0) = x 2 + (T c x + -y (14) 

Defining the critical value a c = f ac ( 0,0) = cr 2 /2 as the 
minimum of this function, there is no solution of eq. (13) 


at a < a c (z.e., no intersection of f Gc with a in fig. 2) 
and hence Dhl(<^ < d c ) = 0. On the contrary, at a > a c 
eq. (13) always has a solution: 


\ADhl t 


(v / 4 (a ff 2 ac) +1 - 1) (a > a c ) 

(15) 

<77 + 0((a- a c ) 2 ) (a > a c ) 


Physically, a c is thus the lower threshold for the coupling 
of the diffusion coefficient with the overstressed regions 
in eq. (11), such that a self-sustainable plastic diffusion 
can be reached even in the absence of an external shear 
rate. However, because of the dissipative nature of plastic 
processes, such a stationary regime cannot be realized 
without some energy input to the system, and it is thus 
an artifact of the model for 7 = 0. In the next section we 
will study the situation where such an external energy 
input will be provided via the macroscopically applied 
steady shear. 


D. Stationary solution at low shear rate 


In the presence of a finite constant shear rate, the sta¬ 
tionary PDF is asymmetric and predicts a finite macro¬ 
scopic stress gm{ 7, <a)- We have at 7 > 0: 


/<r c (x, y) = x 2 + 


+ (\A + + £jj) tanh (^) 

tanh (77 + + 

(16) 


with x = VDt and y = Go'yr/x 2 and, as illustrated in 
fig. 2 for three values of 7 > 0, there is always a solu¬ 
tion of eq. (13) for the stationary diffusion coefficient 
Drl, which can be determined numerically by combin¬ 
ing eqs. (16) and (13). Note that, dimensionally, we have 


[x] = [1 /y\ = [<J C ] and 



= [«] = \ Dt ] = [7] • 


In the limit of vanishing shear rate 7 —0, the lowest- 
order scalings of Dhl(t) can be written down explicitly 
(see Appendix B). Distinguishing the two regimes a ^ a c 
and their limiting case a = d c , one thus obtains: 


a > a c : Dh L t ps D H l( 7 = 0)r (see eq. (15)) 
< a = a c : D nh T « (|j) ' (G 0 7 t ) 4/5 


a < a c : Dn^r « C Gq 7 t 


(17) 


where G satisfies the following implicit equation 

C<7 c tanh (^) = a (18) 

In addition, we get the specific perturbative expansions 






















7 


for a close to or much smaller than a c \ 


a > a c : 

-Dhlt 

(~) - A c ^ 2 


a<a c : 

Dwlt 

1 cr 2 . 

~ /— n= - G07T 

2V6 \J OL c — OL 

(19) 

a « a c : 

-Dhlt 

« — G07T 

CTc 



Note that the perturbative expansion eq. (19) involves 
two limits that cannot be exchanged: the limit 7 —0 
should be taken first, before the limit a a c . Consider¬ 
ing the limit a a c at fixed 7 > 0 in eq. (19) leads to in¬ 
consistencies, as -Drl would then diverge when a —>• a~. 
This divergence is non-physical, since T st r = Dl ^ r < 1 
is bounded by 1 as a probability, and it actually appears 
as soon as the assumption of a ‘low’ 7 breaks down, re¬ 
stricting the validity range of the previous perturbative 
expansion of T>hl(<^ < ol c ) to shear rates at least lower 

than 7 *(<a) = This upper bound goes to 0 in 

the limit a —>• a~. As for the third limit a a c , we re¬ 
cover as expected by self-consistency that, if we remove 
the coupling between the diffusion and the plastic activ¬ 
ity, by setting the coupling parameter a close to zero, the 
diffusion coefficient should vanish as well. 


Once the diffusion coefficient T>hl is known, the 
stationary distribution T s t(cr) is fully determined, and 
its average value gm = da a V s t (&) can be com¬ 
puted, yielding eventually an analytical prediction for the 
macroscopic stress in the limit 7 — 0 : 


a > a c : 


a = a c : 
a < a c : 


O'M 


O'M 


1 + 


i+X'i 


1 


® x o 24 ^ 2 ; f ac ( Xo , 0 ). 
T (G 0 7t ) 1/5 


Gojr 


2 4 / 5 x 3 3 / 5 
cr M ~ cry(c7 c ) + A(a c ) • (Gq7t) 1/2 


( 20 ) 


with #0 = v/£>hl(7 = 0 )t and / (Tc respectively given by 
eqs. (15) and (14). In the last case < d c , the prediction 
of a threshold stress and of a power-law dependence on 
the shear rate, corresponds to a Herschel-Bulkley law of 
exponent 1/2. The corresponding macroscopic stress is 
given by: 


oy(o c ) = C 


M 2 


Cg c tanh (|^) 
and the prefactor of the power law by: 


A(cr c )J^ coth(g) 


+ 


C 3/2 


1 + 


cosh (j$) — 1 
1 “£sinh(^) 


( 21 ) 


( 22 ) 


with the factor G determined by eq. (17) (see Ap¬ 
pendix B). We can finally give the perturbative expan¬ 
sions for the macroscopic stress corresponding to the lim¬ 
its given in eq. (19), first in the limit of small xo, then 
expanding the hyperbolic tangent depending on whether 
G diverges or tends to zero: 


^ ~ OM ( 14 ) o 2 c (15) oy _ \ —2 

a z, a c : ^ ~ ^ - (<a - a c ) 

Goyr 12 xq 12 


^ - - ( 21 ) V 2 C 1 /- \ 1/2 

«s«« : ^ = 12c “ yg ( “'- o) 


A 


(22) C 3 / 2 (19) 


2 3/2 x g3/4 (“<= 


-3/4 


_ (21) <7 C (19) (J c a 

a<Sa c : cry « — - C « —- 

2 2 (J ^ 


(22) y/C (19) 1 / a \ 1/2 


2 Va, 


(23) 


Note that although gy tends to zero in the limit a a~ 
and is thus physically well-behaved (predicting the dis¬ 
appearance of the macroscopic yield stress close to d c ), 
the validity range of its expression has been fixed previ¬ 
ously in the perturbative expansion of the diffusion co¬ 
efficient, with the upper bound 7 * (a) ~ (a c ~ a) 1/2 . A 
similar validity range can be defined in the limit a , 

when Xq becomes comparable to Goyr, with another up¬ 
per bound 7 **(<a) = 12 ^/~ 4 ^ squeezing this regime in 
the vicinity of a c . As for the third limit a < a c , we find 
that if we progressively remove the diffusion term (z.e., 
£ p i (t) = 0 in eq. (5) of our toy model), we eventually de¬ 
stroy the Herschel-Bulkley behavior and the macroscopic 
yield stress tends to g c / 2 , which is the arithmetic average 
between the local yield stress g c and the local stress after 
a plastic event (0 in the full relaxation assumption). 


One of the strengths of the HL model is that it can 
predict three different scaling regimes for the station¬ 
ary macroscopic stress gm{ 7 ) ^ 7^ at low shear rate, de¬ 
pending on the intensity of the plastic events feedback on 
the distribution of local stress, a feedback quantified by 
a = D UL /T st . Note that, on the one hand, a complete 
study of the existence and uniqueness of the stationary 
solution in the limit 7 —0 can be found in refs. [23-25], 
in a dimensionless formulation. On the other hand, the 
only three scalings given explicitly in the original paper 
by Hebraud and Lequeux [18] are partial expressions for 
gm{ 7 ) at a > d c , a < a c and a = a c . However, the ex¬ 
plicit dependence of all the predictions with respect to 
the fixed value of threshold stress g c is crucial for the 
disordered generalization of the model in sect. IV D. 


























E. Connection with the KEP model 


The original HL model [18] defined by eqs. (8)-(ll) 
thus proposes a simplified mean-field mesoscopic scenario 
that can reproduce different rheological laws for cr M ( 7 ) 
depending on the value of a. However, it does not pro¬ 
vide any justification for either the diffusion contribution 
Dhl(*) in the evolution equation (8), or for the lin¬ 
ear closure relation (11) relating the diffusion coefficient 
to the plastic activity with the control parameter a > 0. 
The KEP model [27] precisely provides such a justifica¬ 
tion, and as such we briefly recall its assumptions in order 
to complete the framework presented in this section. 

Providing a derivation of the HL model from a spatial 
mesoscopic picture, the KEP model actually allows for 
the description of spatial inhomogeneities and the spa¬ 
tial propagation of stress after a plastic rearrangement. 
It thus starts from an evolution equation d t Vi(cr,t) sim¬ 
ilar to eqs. (8)-(9), but instead of the diffusion term, it 
includes an explicit spatial (discrete) dependence i and 
a nonlocal operator (denoted ‘£('P,'P)’) coupling distant 
regions through an elastic propagator. Note however that 
the local distributions Vi(cr,t) are assumed to be inde¬ 
pendent in this operator, allowing for the factorization 
of their joint distributions (z.e., V%j = ViVj). 

The KEP model eventually transforms into an effective 
HL equation, with an effective local shear rate 7 i(t) and 
a diffusion coefficient Di(t ), using the following set of as¬ 
sumptions: (i) A plastic rearrangement occurs at a site j 
soon after its local stress a' has exceeded the local yield 
stress, so |cr'| ~ cr c ; (ii) The stress is then fully relaxed at 
site j ( a' —y 0 ), so that the stress fluctuation that prop¬ 
agates from site j to site i is Sa^ « —G^-cr' « —G^cr c , 
with Gij the microscopic elastic propagator; (in) This 
stress 5(j\^ acts as a perturbation that can trigger a plas¬ 
tic event only if the site i is already on the verge of yield- 
§a U) 

ing, z.e., |-^—| « \Gij\ <C 1. In other words, assuming 
either that a plastic rearrangement occurs sufficiently far 
away, or that the amplitude of the elastic propagator is 
small enough, then it can be treated as a noise acting 
on the mean-field local stress, in the spirit of eq. (5)-( 6 ). 
The KEP model then relates the effective diffusion coeffi¬ 
cient D to the plastic activity T, both at a given position 
r, according to: 


D( r, t) = rh(a c ) dlT(v, t) + a(a c ) I\r, t) 
j rn{(7 c ) = b 2 a 2 G 2 nn 
\ & ( a c) = 7 E<(#) Gij 


(24) 

(25) 


with b the discrete lattice parameter, and G nn the 
nearest-neighbor propagator, that is the elastic propa¬ 
gator between neighboring blocks. In the homogeneous 
and stationary case, the linear closure relation ( 11 ) is re¬ 
covered, and the coupling parameter a can be identified 
with a(a c ) given in eq. (25). So the diffusion term in the 
HL model is justified and properly related to the physical 
quantities of the systems, namely the typical value of the 


yield stress cr c and the microscopic propagator Gij. 

The structural disorder is not included explicitly in 
the HL or KEP models, since the local potential energy 
landscape - denoted Vq(J) in sect. II - is summarized 
into a single value of the yield stress cr c . In the next 
section, we will consider the generalization both of these 
constructions, by including a distribution of yield stress 
values p(cr c ). 


IV. HEBRAUD-LEQUEUX MODEL WITH 
STRUCTURAL DISORDER 


The assumption of a single value for the yield stress 
a c is of course restrictive, and a natural generalization of 
the HL and KEP models relies on the existence of an a 
priori distribution of such threshold stresses p(cr c ). In the 
physical picture presented in sect. II, it can formally be 
derived from the distribution of the local potential energy 
landscape Vo(^), more specifically from the distribution 
of its inflection points, as illustrated in fig. 1 (right inset). 
Our motivation for studying a generalization of the HL 
model that includes a structural disorder explicitly, via a 
distribution of yield stresses, is in particular to examine 
the robustness of the HL predictions for the macroscopic 
stress at vanishing constant shear rate 7 —>• 0 : 

{ a > a c : D ~ 7 0 => cr M ~ 7 

a = ol c : D ^ 7 4 / 5 => o M ~ 7 1 ^ 5 (26) 

a < a c : D ^ 7 => a m = cry + A 7 1 / 2 

The questions we wish to address are the following. Do 

we still predict three regimes controlled by the parameter 
a, with these specific exponents? If we predict qualita¬ 
tively the same behaviors, what are the corresponding 
prefactors? How to define the ‘critical’ value a c l And 
can we still derive such a disordered HL model from a 
KEP-like construction? 


A. Definition of the disordered HL model 


Our disordered HL model is defined by the following 
evolution equation for the joint PDF V(a c , cr, t) of having 
a local stress cr and a local yield stress cr c at time £, under 
an external shear rate 7 (£): 

d t T(<j c , <r,t)=- G 0 j(t) d a V + D(a c , t) d^V 

- i / hl(ct, <7 c ) V + T(t) 5(a) p(a c ) 
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with the following specific rate z^hl and the plastic activ¬ 
ities T(a c ,t) and T(t): 

^hl(ct,(T c ) ( = -6{\a\ - a c ) (28) 

T 

poo ~ 

T(t) = (^hl(ct, cr c )) f{ac ^ t) = J da c r(cr c , t) (29) 

f(a c ,t) = - [ da'V(a c ,a',t) (30) 

T J \<j' | >cr c 

where 0 is the Heaviside function and S the Dirac distri¬ 
bution. The evolution equation (27) differs from eq. ( 8 ) 
only in its last term, which includes the a priori distri¬ 
bution p(cr c ). It corresponds to the hybrid stochastic 
dynamics defined by eq. (5)-(4) for the mean-field local 
stress, but randomly selecting according to p(a c ) a new 
yield stress value cr c , after each plastic stress release in 
eq. ( 6 ) (in the same spirit as in the SGR model [14, 15] 
and its precursor the Bouchaud trap model [28]). The 
choice (28) assumes that the rate of plastic events can 
be approximated by a fixed value 1 /r in a locally over¬ 
stressed region (where \a\ > cr c ). 

We can now distinguish the global plastic activity T(t) 
from its partial counterpart T(cr c ,t), defined respectively 
as the mean rate of plastic events for all the system, and 
the rate restricted to a given value of yield stress. Sim¬ 
ilarly, we can distinguish the joint PDF V(cr c ,cr,t) from 
the PDFs of local stress and local yield stress, respec¬ 
tively: 


V(a,t) = 

poo 

/ da c V(a c ,a,t ) 

Jo 

p OO 

(31) 

P( a c,t) = 

/ daV(a c ,a,t ) 

J —OO 

(32) 


The normalization J^° oo daV(a,t) = 1 imposes the first 
equality in eq. (29). By integrating eq. (27) over a un¬ 
der the assumption that both V(a c , cr, t) and d a V(a c ^ cr , t) 
vanish when |cr| —>> oc, we obtain the additional relation, 

d t p(cr c , t) = T (t) p(a c ) - f (cr c , t) (33) 

This equation describes how the distribution of local yield 
stress evolves in the system, depending on how the sam¬ 
ple has been prepared (the initial condition p(cr c , 0 )) and 
how fast plastic events refresh the potential energy land¬ 
scape. This picture simplifies only in the stationary state 
at fixed shear rate, where 

dtp s t(cr c ) = 0 => f st (a c ) = r st p(a c ) (34) 

The partial activity T st (a c ) is then proportional to the 
a priori distribution p(cr c ) (which in general differs from 
the dynamical steady-state distribution p st (cr c )). 

As for the diffusion coefficient, it can generically be 
denoted D(a c ,t), allowing for a different diffusion coeffi¬ 
cient for each value of a c . Such a dependence on a c would 
mean that the mechanical diffusion stemming from the 


plastic events would be different depending on the yield 
barrier to overcome, z.e., modifying the proportion of ac¬ 
tive sites depending on the barrier height. This case can¬ 
not be completely ruled out physically, and it will be 
briefly discussed in Appendix C. Nevertheless, a gener¬ 
alization of the KEP construction - that we will present 
at the end of this section - rather suggests that the me¬ 
chanical diffusion coefficient should be assumed to be the 
same for all regions of the system, disregarding the local 
yield stress value, in the sense that it should gather in one 
parameter the collective feedback of all the overstressed 
regions that might yield. In that case, since a plastic 
rearrangement with a higher barrier will release a larger 
stress in the rest of the system, we will assume that the 
diffusion coefficient D(a c ,t) will simply be replaced by 
the generalized linear closure relation: 

poo ~ 

£>hl(£) = / da c a(a c )T(a c ,t) (35) 

Jo 

where a(a c ) is an ad hoc set of parameters of the model. 
In the stationary case, we recover the same closure rela¬ 
tion (11) as in the original HL model, by using eq. (34): 

^HL = ^effr st (36) 

poo 

<a e ff= / da c a(cr c ) p(a c ) = (a(a c )) (37) 

Jo 

Note the introduction of the notation (0(a c )) for the dis¬ 
order average over p(cr c ) of an arbitrary observable 0(a c ). 

The expressions for our disordered HL model will 
be given for a generic distribution p(cr c ), and the per¬ 
turbative expansions for the diffusion coefficient and 
the macroscopic stress will be valid as long as the 
moments of this distributions are finite (which is al¬ 
ways the case in physical systems). In particular, the 
choice ps(crc) = d(a c — crj) corresponds to the standard 
HL model. However, whenever needed for explicit com¬ 
putations regarding the graphs, we will consider specifi¬ 
cally an exponential distribution for the threshold energy 
E — similarly to the structural disorder included in 
the SGR model [14, 15]. In terms of the distribution of 
<r c , this corresponds to 

po(a c ) = 2 a c exp (-<7^) (38) 

with the mean value (a c ) = y/n/2 and the second mo¬ 
ment {&c) = 1- Rs higher moments can be computed 
straightforwardly, and their rescalings with respect to 
the mean value / ( cr c ) k yield only constant factors 
of order 1. 


B. Stationary solution at fixed shear rate 

We focus exclusively on the stationary solution at con¬ 
stant shear rate 7 , as it is the generalization of the HL 
predictions recalled in sect. III. Thanks to eq. (34) the 
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partial plastic activity can be straightforwardly related 
to its global counterpart with the distribution p(cr c ). As¬ 
suming that the stationary PDF P st (( 7 c,(j) exists, the 
corresponding global plastic activity r st and stationary 
diffusion coefficient D are well-defined. Hence the deter¬ 
mination of the stationary solution of our disordered HL 
model proceeds in the same way as in the original HL 
model (see sect. III). The corresponding explicit expres¬ 
sions are given in Appendix D, along with some technical 
hints regarding their derivation. 


For a given cr c , the stationary joint PDF V s t( cr c cr ) 
can be written as P st (cr c ,cr) = Pcr c (a) p(a c ), which de¬ 
fines V ac (a). Using eq. (27) in the stationary case, we 
find that the distribution Verier) obeys, for a fixed <r c , 
an equation of the same form as that of the standard 
HL model, namely eq. ( 8 ). However, the distribution 
Vcr c ((j) is not normalized to 1 , but instead it satisfies, 
using eq. (32) 


r^>)=^ ( 39 ) 

J— oo PKPc) 

Then, using the same solution procedure as in Sect. IIIB, 
we end up with the relation 


Lc 

Psti^c) = r st rp(cr c )- 


(y/Dr 

yv-LJT, Dr j 


(40) 


where f ac is exactly the same function as in eq. ( 12 ), for 
instance the parabola (14) in the absence of shear rate 
and the function (16) in the presence of a constant shear 
rate. We emphasize the key role that will be played by 
this function f Gc in the present study of our disordered 
HL model. Integrating this last expression over cr c , we 
obtain the counterpart of eq. ( 12 ) for the global plastic 
activity, 


r A (VDt, g $ i ) 

r st r J da c p(a c )- V — -- = 1 . ( 41 ) 

Since we have restricted ourselves to the case where the 
stationary diffusion coefficient takes a fixed value D in¬ 
dependent of < 7 C , as in eq. (35), the previous relation sim¬ 
plifies to 


r stT ■ 


/eff 


tj Go7r\ 

Dr J 


Dr 


= 1 


defining the following effective function, 
f eS (x,y) = (f„ c (x,y )) 


(42) 


(43) 


Eq. (41) can be used to compute r st , at least numerically 
if not analytically, for any choice of p(cr c ) and D. Com¬ 
bined with the closure relation (36), it provides us with 


the generalized counterpart of eq. (13): 

U (vs;, q ^L) = «.» (44) 

from which D = Dhl can be determined uniquely as a 
function of the shear rate 7 and of the effective coupling 
parameter a e ff. 

We emphasize that all this procedure is again quite 
generic - with respect to the choice of the rate z/(cr, <j c ) 
and of the closure relation - and that it has the same 
geometrical interpretation as the one illustrated in fig. 2 . 
The more generic case of a diffusion coefficient that would 
depend on cr c , as initially included in eq. (27) for the 
evolution d t V(cr c , cr, £), is discussed in Appendix C. 

Moreover, in the stationary case and with a diffusion 
coefficient independent of cr c , we have direct access to 
the distribution of local yield stress values, by combining 
eqs. (34), (40) and (43): 


Uc 

P s t(o- c ) = />((J C ) —j~Z 

\fo c 

confirming that the distributions p st (cr c ) and p(cr c ) differ 
in general, except if p(cr c ) reduces to a Dirac distribution. 
This generic expression further simplifies with the closure 
relation (44), replacing the denominator with a e s and 
fixing D H l = D( 7 , a e ff)- 

In fig. 3, we have illustrated the stationary joint distri¬ 
bution V st (cr c ,<j) for different values a c (which have the 
same functional form in the stationary case) and the cor¬ 
responding complete distributions p st (cr c ) and V s t (cr) for 
the specific distribution po(a c ) given in eq. (38). The be¬ 
havior of the distribution p s t (cr c ) depending on a and the 
shear rate 7 will be examined in sect. IV E. But for now 
we will consider the diffusion coefficient and the mean 
stress in the stationary case, following the same struc¬ 
ture as for the standard HL model in sect. IIIC-IIID. 


G 0 jr \ 
Dr ) 


D Gq7t \ 
Dr I 


(45) 


C. Stationary solution in the absence of shear rate 

The unsheared case (7 = 0) can be used as a bench¬ 
mark for the comparison between the disordered HL 
model and its original counterpart. 

The stationary joint PDF is symmetric with respect to 
< 7 , and it thus predicts as expected no macroscopic stress. 
But more importantly, applying the definition of / e fr on 
eq. (14), we have: 

X m '/Dr, /eff (x, 0) = x 2 + (a c ) X + i (<Te) (46) 
ac = /eff(0,0) = i (al) (47) 

where all the factors ~ a k in eq. (14) have been trans- 
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oy, resulting only in slight quantitative modifications of 
the predictions. 


D. Stationary solution at low shear rate 


In the presence of a finite constant shear rate, the sta¬ 
tionary joint PDF is asymmetric at each fixed cr c , and 
it thus predicts a finite macroscopic stress cfm( 7 , ayff). 
From eq. (43), we can compute at 7 > 0 the function f e q 
for an arbitrary distribution p(oy), at least numerically, 
from the known expressions for f Gc defined in eq. (16). 



Figure 3. Top: Stationary joint distribution V(cr c ,cr) 
at a e ff = 0.3 < Oi c and Go'yr = 0.1, for four fixed values 
cr c G {0.5,1,1.5, 2}. Each slice of the distribution has the 
same functional form at each cr c , with an asymmetry due to 
the finite shear rate 7 , and is normalized by definition to 
Pst(cr c )- Bottom left: Corresponding dynamical distribution 
of local yield stress p s t(cr c ), as defined in eq. (32), with the 
underlying a priori distribution po(cr c ) of eq. (38) plotted as a 
blue dashed line. Bottom right: Corresponding distribution of 
stress 'Pst(cr), as defined in eq. (31); its precise shape depends 
on po(a c ). 


formed into moments (cr k ). With the updated definition 
(47) of the ‘critical’ value ay, we recover the same two 
regimes for the diffusion coefficient as before, respectively 
— 0 at Qjgff Qig and, at ojgfj ^ ojg. 


For the disordered HL model (27)-(30), we need the 
following definition of ‘critical’ values for the coupling 
parameter: 

«c(<7 c ) = \ a l > a c = (a c ((Tc)) = I 7c} (49) 

Combining it with the stationary closure relation (36), 
we eventually obtain, in the limit of vanishing shear rate 
7, the following lowest-order scaling of Dhl(t)? as § en " 
eralizations of eqs. (17)-(19)-(20) (see also Appendix B): 

' ayff > 07 : D ul t « D H l(7 = 0)t 
« eff = OL c : D kl t « C (Go7t) 4/5 

< 

ayff < ay : Am t ~ C Go 7 T , 

with (Gcr c tanh (f^)) = ayff. 

(50) 

This last equation, which implicitly determines the pref¬ 
actor C at a e ff < 07 , can be rewritten as: 

t anh fe) _ agl \ = „ (J1) 

2 C a c(^c) j 



(see eq. (48)) 


Cm 


24(<j c ) 


2/5 


7-Dhlt 


h K) 71 - l) (tteff > a c ) 

+ o ((a e ff - a c ) 2 ) (a e ff > « c ) 

(48) 


So we can conclude that, at least in the case 7 = 0 , 
the physical interpretation of ay - as the lower threshold 
for a possible self-sustained plastic diffusion - is robust 
to the addition of structural disorder. Both the func¬ 
tion f Gc (14) and the predicted diffusion coefficient Dhl 
( 15) can actually be straightforwardly generalized by a 
proper averaging over the random values of yield stress 


where we see the competition between terms tanh (u)/u 
and ayff/ay(cr c ), depending on the value of oy. For a 
generic distribution p(a c ) and an arbitrary value of ayff, 
the prefactor C'(oyff) has to be determined numerically. 
The different behaviors of the diffusion coefficient are 
illustrated in fig. 4 (top): at low shear rates it scales 
according to eqs. (50), whereas at large shear rates we 
have Dhl't = a e ^T st r « ayff, the system being essen¬ 
tially over stressed everywhere. 

We can nevertheless obtain analytically the following 
specific perturbative expansions, which are the counter- 
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parts of eqs. (19): 


C^eff ^ OL c • 
<^eff ^ <^c : 


C^eff . 


^HLT 




(48) f a e ff — a c 


(°c) 

7 c) 


24 (a c - a e g) 


1/2 

Gq-Jt 


-Ohl't « Go-jr 

(Vc) 


(52) 


Note that these two last expressions are obtained by ne¬ 
glecting the contribution of some values of g c in the total 
average of eq. (51). As discussed in Appendix B, we can 
define for a given a e ^ a typical value cr* = 2 C. Close to 
a c , C — y oo and we can safely neglect the contributions 
of g c > cr*, whereas when <a e ff is close to zero, C 0 
and we can neglect the contributions of a c < a*. So, in 
the vicinity of a c = (a c ), we recover in particular the 
same expressions as in eq. (19) and there is an increasing 
quantitative correction of the prefactors the further the 
coupling parameter a e s moves away from a c . These ap¬ 
proximations are of course valid only if we are sufficiently 
close to a c or to 0 , but also if there is a ‘reasonable’ cutoff 
in the distribution p(cr c ) for very large or very small cr c , 
respectively, as is expected to be the case physically. 






Once the diffusion coefficient Dhl is known, we can at 
last compute the corresponding macroscopic stress gm 
using the explicit expressions eqs. (D7)-(D9)-(D10) given 
in Appendix D, as illustrated in fig. (4). The overstressed 
regions (|cr| > a c ) contribute only linearly in 7 to the 
macroscopic stress, and hence it is only the asymmetry 
of the PDF of the understressed regions (|cr| < a c ) that 
can modify the dominant scaling in 7, at least in the van¬ 
ishing shear rate limit that we want to consider. We can 
derive analytical predictions in the limit 7 —)► 0 for the 
three regimes in a e ff, starting from the case at a e Q > a c : 


O'M 


(14) 


1 + 


4x 0 7 c) + 7 c) r 
— 7 ~r~ -v 

24x 0 (fcr c 7o,0) / 

4x 0 (ag) + (7) 

24 (x% + x 0 (<r c ) + (<j2) /2) 


(53) 


G 07 T 


Figure 4. (Color online.) Top: Stationary diffusion coeffi¬ 
cient as a function of the shear rate, for increasing values of 
a e ff £ {0.01,0.3,0.49,0.5,0.51,1} (blue to red), as indicated 
by the black arrow. The dashed black curve corresponds to 
aeff = ol c = 0.5. At large shear rates we have Dhlt ~ a e ff. 
Bottom: Corresponding macroscopic stress as a function of the 
shear rate, for the same increasing values of a e ff (blue to red), 
as indicated again by the black arrow. At large shear rates 
gm ~ 7 so that the predicted behavior is Newtonian, whereas 
at low shear rates we recover the three regimes in a e ff that we 
have predicted analytically for gm( 7 , a e ff) in eqs. (53)-(56). 
This graph is reminiscent, as expected, of fig. 2 in the original 
HL paper [18], although there are quantitative discrepancies 
due to the structural disorder included via po(a c ) of eq. (38). 


over the distribution p(o c ) of threshold stress values. In 
fact, taking care of the different averages according to 
eq. (D5) (and detailed furthermore in Appendix D), we 
obtain for the macroscopic yield stress: 


with xq = a/^hl (7 = 0)t given by eq. (48), and in the 
two other regimes in a e ^: 


a e ff = ol c : 

om 

^ Oi c . 

Om 


7c> 


3/5 


/ \2/5 

Wc) (n . \ 1/5 

n-( GqIt) ' 


2 4 /5 x 3 3 / 5 (<j2) 

a Y + A (Go 7 t ) 1/2 


(54) 


As for the predicted Herschel-Bulkley behavior at 
a e Q < ol c , the macroscopic yield stress Gy and the pref¬ 
actor A are not simply given by averaging eqs. ( 21 )-( 22 ) 


Gy = C 


7 c 2 )/ 2 

C (g c tanh (|g?)) 


(55) 


with C defined by eq. (51). As for the prefactor A , it is 
more subtly obtained by a Taylor expansion of D hl to the 
next order (whose derivation is sketched in Appendix B) . 
Its complete expression, rather cumbersome, is given in 
eq. (E 2 ) and it simplifies in the usual limiting cases of 
interest a e ^ < a c and a e ^ < a c , as presented thereafter. 
The complete behavior of Gy and A as functions of a e Q 
are illustrated in fig. 5 (see Appendix E for the complete 
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expressions). 

We can finally give the perturbative expansions for the 
macroscopic stress corresponding to the limits given in 
eq. (52), first in the limit of small xo, then expanding 
the hyperbolic tangent depending on whether G diverges 
(c^eff ^ ol c , see eq. (E3)) or tends to zero (a e fi <C <a c , see 
eq. (E4)): 


C^eff ^ OL c • 
<^eff ^ O'c : 


^eff < ^' OLc • 


CM (53) <cr^)/(cr2) (48) (cr)) (<7 c ) 2 /(<7 2 ) 


Go7t 12a;g 12 (o; eff - a c ) 

(a c - a e ff) 1/2 (o'c) ^ 


(Jy 

4 r 
cy 

4 R 


V 6 (a c 2 ) 

(«c - a e ff )“ 3/4 (o"c > 3/4 (^c) 


2 3 / 2 x 6 3 / 4 
Q( c C^eff 


(c?> 


(^c) 

1 _ 7 c) \ / Q eff \ 1/2 

v 2 (<7 C ) 2 i V(°c) / 


(56) 



&M ~ a Y 



Illustrated in fig. 5, these behaviors are qualitatively sim¬ 
ilar to the predictions of the HL model without disorder. 
They actually differ only quantitatively from their coun¬ 
terpart expressions (23) since they involve the moments 
of the distribution p(cr c ) instead of simple powers of a c . 
Consequently, we can transpose the physical discussion of 
their counterpart expressions (23) to our disordered HL 
model, taking into account the quantitative corrections 
due to the distribution p{o c ): (i) ay tends to zero in 
the limit a e ff —)► a~ and is thus physically well-behaved 
(predicting the disappareance of macroscopic yield stress 
and hence of the Herschel-Bulkley behavior at low shear 
rate); (ii) nevertheless, the validity range of its expres¬ 
sion has been fixed previously in the perturbative expan¬ 
sion of the diffusion coefficient, with the upper bound 
7*(<^eff) ~ (OLc - ^eff) 1/2 ; (Hi) finally, a similar validity 
range can be defined in the limit a a+, when x\ be¬ 
comes comparable to G 07 T, with another upper bound 

. / \ 12(a c —a e ff) (<r c ) . 

7 **(<a e ff) = —y — / 4 \ / \ 2 that squeezes this regime m 

Oo T \(J c ) \(Jc) 

the vicinity of a c . As for the third limit o;Ca c , we re¬ 
cover that if we remove the diffusion term, we eventually 
destroy the Herschel-Bulkley behavior, and the macro¬ 
scopic yield stress tends to a c / (a c ), which is slightly dif¬ 
ferent from the arithmetic average between the mean lo¬ 
cal yield stress (a c ) and the local stress after a plastic 
event (0 in the full relaxation assumption). 


However, we can comment furthermore on the valid¬ 
ity range in 7 of the Herschel-Bulkley behavior, and 
hence on the definition of the prefactor A. The con¬ 
sistency of the perturbative expansions of the diffu¬ 
sion coefficient (50) at a fixed low 7 requires that 
AnX^eff < a c ) < AnX^eff = OL c ), which implies in gen- 


Figure 5. Top: Macroscopic yield stress ay and prefactor 
A of the Herschel-Bulkley behavior predicted at a e ff < ^c, 
in the limit of 7 -» 0 (for po(a c ) given by eq. (38)). Their 
behavior is in agreement with the predictions in eq. (56). On 
the left, Gy is evaluated by computing the average stress <tm 
at a very low shear rate G 07 T as 10 -5 (dots), showing a good 
agreement with the theoretical prediction eq. (55) obtained 
by using the diffusion coefficient computed numerically 
(continuous line). On the right, the prefactor A is evaluated 
similarly, using again the diffusion coefficient computed 
numerically and the corresponding prediction for ay. 
Bottom: Test of the validity range of the Herschel-Bulkley 
behavior <jm ~ cry + A (Gojr) 1 ^ 2 , for increasing values of 
aeff e {0.01,0.05, 0.1, 0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.49} 
(blue to red, as indicated by the black arrow) and a e ff = 0.5 
(dashed black). The plateaux at low shear rate define the 
prefactor A, and their validity range actually shrinks when 
a e f£ tends to <a c , in agreement with eq. (57). 


eral that G 07 T < (G/G) 5 , and when a e ff < a c that: 


G 07 T < 



( 24 \ 1/2 ({a c -a eS f / 2 \ 

\K)J \ (v c ) 2 ) 


(57) 


This condition defines a much more restrictive upper 
bound for 7 than the condition Arn'r < <a e ff, and it ex¬ 
plains in fig. 5 ( bottom ) the shrinking of the plateaux 
when <a e ff —)► a ~. Given the fact that the scalings of the 
diffusion coefficient with respect to the shear rate are the 
same with or without disorder, this argument is valid 
both for the standard HL model and for our disordered 
version of it. 

So we have shown in this whole section that, when 
we include a distribution of threshold stress values p(cr c ) 
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in the HL model, the three scaling regimes at low con¬ 
stant shear rates for the macroscopic stress gm{ 7 ) are 
the same as in eq. (26). Provided that we general¬ 
ize the definition of the ‘critical’ coupling parameter 
a c = (a c (cr c )) = (cr?) /2, the expressions of the prefac¬ 
tors are almost identical to eq. (23) in the vicinity of 
<a c , the terms being replaced by combinations of the 
first moments of the distribution p(cr c ). We emphasize 
that it is the scaling behavior of the averaged stress 
gm = dcr a V s t{&) which is qualitatively robust with 

respect to a structural disorder; the fluctuations of the 
stress characterized by the full PDF V s t (cr) are of course 
modified, as illustrated in fig. 3. 


E. Stationary distribution of local yield stress 

Pst (cr c ) 


Combining eqs. (45) and (44), we have access to the 
stationary distribution of local yield stress: 


Pst(crc) = p(°c) 


L c 


( WDr 

<a e ff 


(58) 


with f Gc given by eq. (16) and D = Dhl( 7, a e ff) deter¬ 
mined by solving eq. (44). In what follows, we denote by 
O the average of an observable O over the distribution 
p st (cr c ), while (O) still denotes the average of O over the 
distribution p(cr c ). The main features of this distribu¬ 
tion can then be characterized by its mean value and its 
second moment, respectively: 

00 

da c a c p st (a c ) (59) 

OO 

da c alp st (a c ) (60) 

We emphasize again that the stationary distribution 
p st (cr c ) is not equivalent to the a priori distribution 
p(cr c ), as we can show explicitly in the limit of van¬ 
ishing shear rate (see the last remark in Appendix B). 
The set of predictions given in this section are further¬ 
more illustrated in figs. 6 and 7 for the specific choice of 
PoLc) = 2<t c exp (-of). 



cceff ;$ ol c = (of) /2: 


lim 

7—>-0 


Pst (cr c ) 

PLc) 


lim (j c « 
7^0 


lim a^ ~ 

7-^0 


7 f 1 «efA oj 

2a c \ a c ) (<r|) 

77 (, 2a eff \ 77 

7c) V LDJ Li) 
Li) (, 2a eff \ 77 

ld v ld ; {<) 


and secondly the case at a e ff <a c : 

Pst(cT C ) Oc^°) CTc a c ( a eff<«c) G c 


lim 

7^0 p(cr c ) 


lim <j c 

7—>-0 


2a e ff cr, 

77 


7c) 


— 7<=) + 


Lc) 

Var(cr e ) 

7c) 


lim cr;? 

7—>-0 



(62) 

(63) 

(64) 

(65) 

( 66 ) 

(67) 


with Var(cr c ) m (cr;?) — (cr c ) 2 . eq. (66) shows that, when 
the diffusion is suppressed (a 0), we have a~ c ^ (a c ) 
and the distribution p st (cr c ) does not coincide with p(cr c ), 
except if p(cr c ) has a zero variance, in other words if cr c 
can take only one value. 

We consider now the case at a e ff = ct c , for which we 
can write down the two lowest orders at low shear rate: 


Pst 7c) (77°) _7_ J_ 
p(cr c ) 2 a c a c 



with C = 


ML 

24<a c ) 


2/5 

the prefactor of the diffusion co¬ 


efficient, as defined in eq. (50), leading to: 


_( 7 —> 0 ) (of) 

CTc “ LI) 



7c 5 ) 1 ) 

12 LI) C 2 ) 


(G 0 7t) 2/5 (69) 


-2 M 0) 77 

c ~ Lc) 


77 71/2 

(7) ' 


7c 6 ) U 

12 7c 2 ) C 2 J 


(Go7r) 2/5 

(70) 


We consider at last the case at a e s > a c , for which we 
can also write down the two lowest orders at low shear 
rate: 


Art 7c) (7^0) 7o + <JcX 0 + qf/2) + 0 /^ 2 \ , ?1 * 

p(c r c) cr e ff 


We start with the case a e s < <a c , whose strictly van¬ 
ishing shear rate limit predicts: 


with xo = \JDyllt given by eq. (48), Dhl being the dif¬ 
fusion coefficient at zero shear rate. This leads to: 


lim 

7—>-0 


Pst {°c) 

pLc ) 


a 2 tanh (<t c /<t*) 
2a e ff (j c /cr* 


< = 2 G (61) 


_ (M» fe)»S + W>^ + W>/2 + g m (72) 

C^eff 


with C the prefactor of the diffusion coefficient de¬ 
termined by eq. (51). Using the relations (52), we 
can write explicitly the usual specific limits, first at 


^2 7c 2 ) ^o + 77^q +7c)/ 2 + 0 ^ (73) 

^eff 
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In the case a e s > a c , the expansion (71) is actually valid 

only for shear rates such that Gn'yr < 12 ^* c ( a c) ag 
J U/ GqT (<T’) \<j c ) 


already noticed after eq. (56) for the mean stress gm 
expansion above a c . Furthermore, reading directly the 
limit 7 —)> 0 of eqs. (71)-(73), we can write explicitly the 
expansions at <a e ff > ol c using that xq « Qe r? c : 


lim 

7—>-0 


Pst (CT C ) 

pWc) 


lim g c 

7^0 


lim g% 

7^0 


G c ^ G c Ct e ff G c 
2 Qy (tT^>) oy 

M) + M) _ d 

<*t> + (<7 C > V «c ; 

( ff c) 7c) [aeB 7 

(°c) (<^c) V «c ) 


(74) 

(75) 

(76) 


g c Var p(t 7 c ) 



Figure 7. (Color online) Mean and variance of the dy¬ 
namical distribution of local yield stress p(cr c ), as a func¬ 
tion of G 07 T, for three values of o e ff- a e ff = 0.3 (dotted 
blue), a e ff = 0.5 (continuous black) and a e ff = 1 (dot-dashed 
red). The corresponding curves are sections of fig. 6 at fixed 
a e ff- In the limit of large shear rates, we recover as ex¬ 
pected from eq. (77), for po(a c ), that oy ~ (cr c ) 0.886 and 
Varp(cr c ) ~ Var p (cr c ) ~ 0.215. 


1 



Figure 6 . (Color online) Mean and variance 

of the dynamical distribution of local yield stress 
p(a c ) as a function of a e ff, for increasing values of 
G 0 7 t e {lO -7 ,10 -6 ,10 -5 ,10 -4 ,10 -3 ,10 —2 ,10 -1 } (blue to 
red), as indicated by the black arrow. The mean and variance 
are computed from eq. (58). The expected limits o e ff 0 
and a e ff o c , computed at 7 0 for po(a c ) of eq. (38), are 

recovered. Results for larger shear rates are shown in fig. 7, 
for fixed values of a. 


The complete behaviors of the mean value o~ c and the 
variance Var^oy) = cq? — dy 2 are illustrated in figs. 6 
and 7, respectively at fixed 7 and at fixed ayff, for the ex¬ 
ponentially decaying distribution po(cr c ) of eq. (38). We 
can immediately see that a e s = ay plays a special role in 
the limit 7 —» 0 , as it corresponds to the maximum value 
of dy. Physically, the dynamical distribution p st (cr c ) is 
an interplay between the a priori distribution p(cr c ) and 
how fast the local yield stress values are refreshed by the 
plastic events triggered by the shear rate. The larger the 
shear rate, the faster the local yield stress is refreshed 
and the more small energy barriers oc cr 2 can be present 
in the system. On the contrary, the lower the shear rate, 
the less plastic events we have, so the largest barriers 
survive longer and the distribution p st (cr c ) allows larger 
values of a c on average. 

These different sets of predictions for the stationary 
distribution p(cr c ) provide new features to test numeri¬ 
cally or experimentally, in addition to the characteriza¬ 
tion of the mean stress gm{ 7 ) and more generally of the 
complete stress distribution V st (a). The complete map¬ 
ping between the dynamical distribution p(oy) and the 
a priori distribution p(cr c ) is provided by eq. (58). In 
the case of a low shear rate, the strict limit 7 —)• 0 can 
be computed exactly for the different regimes in ayfj; in 
the specific cases of ayfj ay, a e s < a c and ayfj > a c it 
turns out to depend solely on the first moments of the 
distribution p(a c ). Moreover, the next order in Gq'Jt is 
available at ayff = a c and at a e ff > ay, and the discrep¬ 
ancy between these scalings - respectively (Go7t ) 2/5 
and ~ (Go 7 r ) 2 - provides an additional way to probe in 
which regime in oy^ we might be. 

In a numerical or experimental test of these predic¬ 
tions, if we have access simultaneously to the complete 
dynamical distribution p st (cr c ) and to the parameters 
{Go, D, 7 , r, <a e ff}, then we can use eqs. (58), (61), ( 68 ) 
or (71) in order to determine the a priori distribution 
p(g c ). If we do not have access to the complete p st (oy), 
but only to its mean value and variance for instance, then 
the connexion is still possible but we need to assume a 
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given shape for p(a c ) [31]. 

We can also mention the opposite limit of large shear 
rate, in which T st r « 1 (all the sites are brought above 
the local yield stress in the time interval r, in the station¬ 
ary state), so DhlT ~ <^ e ff and we have at lowest order 
in l/y: 


Pst{cr c ) ( 7 ^ 00 ) fcr c ^\A*effj 


Gpjr \ 

C*eff ) 


( 16 ) 


pip c) 


C^eff 


Goyr 


(77) 


Consequently, we predict at large shear rates that 
Pst(< 7 c ) ~ pi^crc)•, with the same mean value and variance, 
as is indeed the case in fig. 7. Physically, at large shear 
rates the system is essentially overstressed and on the 
verge of yielding everywhere, so first the diffusion co¬ 
efficient tends to a constant, secondly the mean stress 
cr M ~ 7 has a Newtonian behavior, and thirdly the dis¬ 
tinction between the dynamical and the a priori distri¬ 
butions of local yield stress a c has been washed out. 

It is interesting at this stage to try to provide a phys¬ 
ical interpretation of the results obtained for p s t(cr c ), or 
in other words, to try to understand how the a priori 
distribution p(cr c ) is reweighted to give p st (cr c ). A naive 
reasoning suggests that if the dynamics of a is dominated 
by the drift 7 , the time needed to go from a = 0 to cr c 
should be proportional to cr c , so that one would expect 
p st (cr c ) oc a C p(cr c ). Similarly, if the dynamics is domi¬ 
nated by the diffusion, the time to go from 0 to a c should 
scale as cr?, so that one expects p st (cr c ) oc g 2 c p(a c ). 

The first situation (dynamics dominated by the drift) 
should occur for large 7 . However, eq. (77) shows that for 
7 —)> 00, p st (cr c ) = p(cr c ), so that the a priori distribution 
is not reweighted by a c as expected from the above naive 
argument. This is actually due to the fact that a does not 
jump to 0 exactly at cr c , but has only a probability 1 /r 
per unit time to jump to 0 when a > a c . For large 7 , the 
local stress a thus reaches values much larger than the 
local yield stress a c before jumping to 0. The ‘lifetime’ 
of a state with a given a c is r + <j c /(Go 7 ), and it leads 
to the reweighting given in eq. (77). 

On the other hand, the diffusive regime is expected to 
occur in the low 7 limit, and for <a e ff > a c so that the dif¬ 
fusion coefficient does not vanish. Here again, the result 
given in eq. (71) seems to rule out the naive expectation, 
since the reweighting of p(cr c ) is not proportional to cr?. 
However, a closer look actually shows that intuition was 
not wrong. Considering a fixed diffusion coefficient D 
(z.e., discarding the closure relation eq. (36)) and taking 
the limit r —>> 0 to make the stress jump to zero sharply 
at cr c , we indeed get from eq. (71) that the reweighting of 
p(cr c ) is proportional to o\. The generic deviation from 
this scaling in eq. (71) again results from the fact that 
!t is finite, so that cr can increase above cr c . This can be 
quantified by comparing r to the diffusion time D to 
reach cr c starting from a = 0 : the small r regime corre¬ 
sponds to r <C cr^/D, or equivalently to Dr <C cr 2 C . Note 
that taking into account the closure relation eq. (36), the 
self-consistent diffusion coefficient Dhl is determined via 


the product T>hiT (see eq. (50)), so that the limit r —)> 0 
cannot be taken at fixed T>hl- Actually, the control pa¬ 
rameter is a e ff rather than r, and one finds that DhlT is 
small for a e s close to a c (with a e s > a c ); in this regime, 
one recovers a reweighting proportional to o 2 c to leading 
order, as seen in eq. (74). 


F. Connection with the disordered KEP model 


Now that we have generalized the analytical predic¬ 
tions of the original HL model by including the distribu¬ 
tion p(cr c ), we can conclude this study by examining the 
disordered extension of the KEP construction [27] (re¬ 
called in sect. HIE), in order to provide a justification 
for our disordered HL model, and in particular for the 
choice of a diffusion coefficient independent of cr c . 

First the evolution equation dtVi(a^t) in [27] can be 
extended to dtVi(a c , cr, £), allowing for stress propagation 
between regions with different yield stress values, replac¬ 
ing for that purpose the operator 4 ,VY by: 


,a c ,aiV,V) = f dc t' c - ( da' Pj (cr' c ,cr',t) 




Pi (<r c , a + 5a. 


U \t) 


~ Pi (<r c ,pt) 

(78) 


Note that this operator implicitly assumes that the 
local distributions Vi(cr,t) are independent, allowing 
for the factorization of their joint distributions (z.e., 
Vij = Vi Vj ). This evolution equation then transforms 
into an effective HL equation, with an effective local 
shear rate 7 i(t) and a diffusion coefficient Di(t), with 
the following slightly modified set of assumptions: (i) A 
plastic rearrangement occurs at a site j soon after its 
local stress a' has exceeded the local yield stress, so 
|cr'| ~ ac J K (ii) The stress is fully relaxed at the site 
j (a' —)> 0 ) hence the stress that propagates from the 
site j to the site i is ~ —G^-cr' « —Gij<j^\ with 

Gij the microscopic elastic propagator, (in) This stress 
acts only as a perturbation that can trigger a plastic 
event only if the site i is already on the verge of yielding, 

z.e., | ~fif | ~ | Gij —fij | 1. So this construction remains 

valid assuming either that a plastic rearrangement occurs 

sufficiently far away, or that the amplitude of the elastic 

O') 

propagator is small, but only as long as the ratio 

cr c ’ 

remains sufficiently small as well, a criterion that might 
constrain the variance of the distribution p(cr c ). 

With the operator Ci^ ac , a (V,V) being a linear combi¬ 
nation of the contributions of the different values of cr c , 
the KEP relation between the diffusion coefficient D( r, t) 
and the local plastic activity T(r, t) of eq. (24)-(25) be- 
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comes: 



rh(a c ) d^T(cr c , r, t) + a(a c ) r(<r c , r, t) 


(79) 

and in the stationary case, using eq. (34), we finally ob¬ 
tain: 


^kep(i*) = m e ff [p\ d^T st (r) + [p] r st (r) ( 80 ) 

/ m eff [p] = / 0 °° da c m(a c ) p(a c ) = (m(a c )) p , , 

\ aeffM = / 0 °° da, c a(a c ) p(cr c ) = (a(a c )) p 

This last result justifies, as anticipated, our choice of the 
time-dependent closure relation (35) in general, the linear 
closure relation (36) and the definition of the coupling 
parameter (37) in the stationary case. 


V. DISCUSSION AND OUTLOOK 


A. Discussion on the HL model assumptions 


Regarding the definition of the original HL model and 
its disordered counterpart, one pending issue is the in¬ 
terpretation of the physical mechanism underlying the 
choice of the plastic rate phl(<T & c ) = y#(|cr| — cr c ) in 
eqs. (9) and (28), i.e ., the assumption of a typical fixed 
rate 1 /r of having a plastic event when the local stress 
exceeds the local threshold a c . It has been highlighted 
in an earlier work, that one other key ingredient in order 
to obtain a shear-rate dependent flow-curve in athermal 
systems is the existence of at least one additional intrin¬ 
sic timescale, that will be a material-dependent prop¬ 
erty [17]. This timescale has been identified as the dissi¬ 
pative time, which describes roughly the typical duration 
of the local relaxation process. The original definition 
of the fixed rate 1/r in the HL model is not equivalent 
to this dissipative time, but it rather introduces a ‘lo¬ 
cal overshoot’ regarding the local yield stress value. The 
physical interpretation of this process remains somehow 
unclear. However, it is possible to interpret the time r in 
the HL model as the dissipation time defined in ref. [17] 
by introducing a small correction in the definition of the 
macroscopic stress, namely 


^.corr — 
a M = 


= da c damm(a,a c )V st (cr c ,a) 
JO J R 

(<r c ) Dt 

Gotl\\ 

’ Dt )/ 


(under) 
= (J M 


+ 


(82) 


in other words assuming that in the ‘overstressed’ re¬ 
gions the local stress does not exceed its local yield value. 
Such a modification of the macroscopic stress definition 
does not alter the rheological behavior of the HL model 
at low shear rate, at least for the regimes a e ff < a c (at 
<a e ff > a c it induces a finite macroscopic yield stress). For 
small driving shear, this average is in fact dominated by 
the contributions of local stresses a below the local yield 


stress values cr c , and the Herschel-Bulkley exponent is 
thus robust with respect to this subtle change. Note how¬ 
ever that the large shear rate regime will be influenced by 
such a correction, so we think it is important to eliminate 
the rather unphysical existence of locally overstressed re¬ 
gions for future considerations. Also the existence of a 
finite dissipation time allows for the appearance of shear 
localization [29, 30], a feature which is however out of the 
scope of the present study. 

Another very strong assumption of the HL model is 
the full relaxation of the local stress after yielding, as in 
eq. (5) of our toy model. This scenario is not always - ac¬ 
tually, rather rarely - satisfied [31]. Partial relaxation is 
expected to modify the prefactors in the predicted scal¬ 
ings, for which we are still missing explicit analytical ex¬ 
pressions in this case. Since in this study we have focused 
on the mechanism that leads to the onset of non-linearity 
and on its corresponding exponent, we are not concerned 
with the effect of partial relaxation of stresses. For a 
quantitative comparison, we would need either to com¬ 
pute numerically the stationary joint PDF V s t (cr c ,cr) at 
fixed 7 replacing in eq. (27) the S(a) with a more gen¬ 
eral distribution A (a). Or one can alternatively study 
numerically the whole set of physical quantities we have 
defined, starting from the hybrid dynamics of eqs. (5)-(4) 
on a set of independent sites with local stress {cq(£)}. 

Finally, we can comment on the assumptions regarding 
the underlying distribution of local yield stress needed for 
our results to hold. In the definition of our disordered HL 
model, we have assumed a generic p(cr c ) whose moments 
are finite. For instance, the perturbation expansion of the 
diffusion coefficient at low shear rate already involves its 
fourth moment (&£)• However, although the mean value 
and the variance are reasonably robust features of a dis¬ 
tribution, in practice its higher moments can be strongly 
sensitive to the system size and to the limited available 
statistics. So, although the power-law expansions of the 
HL mean stress (hence the qualitative behavior of the 
rheological law at low shear rate) are predicted to be 
robust to the addition of structural disorder, their corre¬ 
sponding prefactors depend on a combination of moments 
(&c) that might display a dependence on the system size. 
A direct comparison between atomistic simulations and 
our analytical predictions should thus take into account 
such a possible dependence, and would require a careful 
characterization of the distributions of local yield stress 
p(a c ) and p(a c ). 


B. Summary and outlook 

In this study we have identified on a mean-field level 
the necessary ingredients for the modeling of yield stress 
materials in the case of athermally activated yielding 
events. Within the HL model, which we argue to better 
represent the underlying physical picture in athermally 
driven systems, we have studied analytically the robust¬ 
ness of the predictions with respect to an additional and 
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usually important physical ingredient, namely the disor¬ 
der in the yield energy barriers (or equivalently, in the 
local yield stress). We find that, although a key ingredi¬ 
ent in the SGR model, a distribution of threshold stresses 
does not modify qualitatively the HL predictions at low 
constant shear rate, thus predicting a universal critical 
behavior at the flow transition 7 ^ 0 . 

The generality of the different analytical expressions 
in this paper, distinguishing the specific cases or limits 
taken, allows us not only to recover all known results on 
the HL model but also to go beyond them. It enables us 
on the one hand to estimate numerically all the relevant 
physical quantities for a given distribution p(cr c ), and on 
the other hand to consider alternative closure relations 
for the diffusion coefficient. 

In future work we would like to address further ques¬ 
tions to render the mean field equations even more consis¬ 
tent with the underlying physical dynamics. It is known 
that, close to the flow transitions, complex dynamical 
heterogeneities in form of avalanches of yielding events 
play an important role in the plastic response to shear 
[30, 32, 33]. So an important issue that remains is to 
better understand how the spatio-temporal correlations 
of the yielding events can be captured, within the mean 
field modeling approach, notably within the formulation 
of the effective noise term. We emphasize that in the HL- 
like models, the effective noise is approximated by the 
Gaussian white noise assumption, whose variance (i.e. 
the diffusion coefficient) is coupled to the plastic activity 
in the system. In this work we have studied exclusively 
the stationary case at fixed low shear rate, but this as¬ 
sumption might be questioned even more in the further 
study of transient regimes or for an oscillating shear rate. 

Further we will be interested in combining the present 
picture and model with the notion of thermal noise and 
activation, in order to mimic a sheared material that 
is additionally subject to thermal activation of yielding 
events, in the spirit of our toy model picture. It would 
be very interesting to understand in detail the role of dis¬ 
order in this combined picture, already at the mean-field 
level, and later on in a complete statistical field theory 
of the stress field. 
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Appendix A: Stationary plastic activity at fixed 
diffusion coefficient 


In sect. Ill B we have defined the function f Gc as the ra¬ 
tio between the stationary diffusion coefficient D and the 
corresponding plastic activity r st (D), for the standard 
HL case of a single value for a c . This definition was mo¬ 
tivated by the specific closure relation Dhl = aY st given 
in eq. (11), implying that either D = 0, or D > 0 accord- 

ing to the equality / <Tc (v^Dr, 7 ?) = a - 

In fact, the function f c r c (%i,X 2 ) is completely fixed by 
the functional form of the stationary PDF, if we are able 
to solve the specific equation ( 8 ) d t V(cr,t) = 0. If this is 
the case, we can then define the normalized PDF p±(cr) 
from V s t(a ^ 0) = T±rp±(cr) and f±^° dcrp±(cr) = ±1. 
This leads generically to: 


Uo VTr, 


Gojr\ 
Dr ) 


Dr 


= 1 


P+( 0 ) 


p+(0) +p_( 0 ) 

P-( °) 

P+(0) +P-(0) 


f 

J -a, 
r& c 

JO 


dcrp-(c 


dcrp + (c 


(Al) 


and is equal to l/r st (D), see eq. ( 12 ). 

In sect. IV B, we have furthermore generalized this def¬ 
inition to the case with an arbitrary distribution of lo¬ 
cal yield stress values p(g c ). The normalization of the 
PDF implies respectively eq. (41) for a generic D(a c ), 
and eq. (42) for a diffusion coefficient independent of 
ci c . In both cases we can use the function f Gc previ¬ 
ously determined for a fixed value of cr c , as explained in 
Appendix D. 


Appendix B: Perturbative expansions of the 
stationary diffusion coefficient 

In sect. Ill D and IV D, we give the perturbative expan¬ 
sions of the diffusion coefficient Dhl at low shear rate, 
respectively for the standard HL model in eq. (17) and for 
its disordered counterpart in eq. (50). The perturbative 
expansion of the mean local stress gm is then straightfor¬ 
wardly obtained by substituting into its exact expression 
(D10) the expansion of Dhl, and expanding the resulting 
expression at small Go 7 t. 

In this appendix, we first recall exact mathematical 
results regarding the perturbative expansion of Dhl in 
the standard HL model. The general structures of the 
perturbative expansions of Dhl and gm are discussed 
in Chapter 2 of ref. [23], and are given explicitly by the 
theorem 4.1 of ref. [25]: 

,^ 0) f O(7 0 ) + O(7 1 ) + O(7 2 ) + ... («>a c ) 

Dhl » { 0( 7 4/5 ) + 0( 7 1 ) + ... (a = a c ) 

[ Oij 1 ) + 0(j 3/2 ) + ... (a < a c ) 

(Bl) 
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and its corollary 4.2 [25]: 

f 0(Y) + 0('y 2 ) + ... (a > a c ) 

a M « < 0( 7 1/5 ) + 0(7 2/5 ) + ■ • ■ (a = a c ) (B2) 

[ O(j 0 ) + 0( 7 1/2 ) + ... (a < a c ) 

Note that in these references, the quantities of interest 

are made dimensionless with respect to the single value 
cr c , with the definitions /x = a/cr 2 , <fi = Dr /cr 2 and hence 
r st r = Dn^r/a = 0//x. These perturbative expansions 
are systematically constructed and well-controlled, with¬ 
out any a priori knowledge of the corresponding conver¬ 
gent series, using ‘asymptotic expansions’ of the station¬ 
ary solution of the PDF [25] . Although such a procedure 
can in principle be generalized to the disordered HL case, 
it is not straightforward. 


We have thus taken a shortcut for the disordered HL 
case, in order to obtain the lowest order of the expan¬ 
sion given in eq. (26), using the exponents in eq. (Bl) as 
guides in the standard Taylor expansions. This shortcut 
consists in assuming the following ansatz at low 7 : 

0 (7—^0) r 

X = D kl t PS Cl (Go-yr) 61 
_ Go 7 r ( 7 ^ 0 ) (Go 7 t ) 1-<51 ^3) 

V ~ a : 2 * ~C[ 

with 0 < $1 < 1. In the limit 7 —> 0, we have three pos¬ 
sible cases: 

$i = l: x —y 0 , y —y 1 / C\ 

0 < $1 < 1 : y —y 0 (B4) 

$i= 0 : x —y x$ , y —y 0 

for which we expand the function f ac (x, y) given in 
eq. (16) at low shear rate. If $1 = 1, we have: 

lim h c {x,l/Ci) = GiCTctanh ( ) (B5) 

If $1 = 0, we have: 

2 

U c po, Goyr/xl) =x% + a c x 0 + y + O ( 7 s ) (B 6 ) 

And if 0 < $1 < 1 , we have: 

Uc (x,y) = y +a c x-^y 2 +x 2 +0 ( xy 2 )+0 ( y 4 ) (B7) 

We can identify which value of $1 is associated to each 
regime of a e s and determine the corresponding prefactor 
Ci , by solving at lowest order the equation deduced from 
the closure relation (35): 

(Jcjc (x, Gq-jt/x 2 )^ - a eS = 0 


We start from the case $1 = 0, that yields: 

Xo + ( a c) ^0 + ^ 7c) = a eff (B8) 

which admits a positive solution of xq = V^hiT only if 
a e ff > a c = \ (cr 2 ), as gi ven by eq. (48). We then turn 
to the case 0 < $i < 1 , where x 2 ~ 7^1 is small compared 
to x ~ /A/ 2 , so we can cancel the two lowest orders pro¬ 
vided that: 

O^eff = ^ ( a c) = a c 

/ 4 \ (B9) 

(a c ) C 1 / 2 (Goyr)^ 2 - ^ (G 0 ir) 2 ~ 2 ^ = 0 

implying that $1 =4/5 and C\ == C as given by eq. (50). 
The last case $1 = 1 should thus correspond to a e ff < a c, 
and it actually yields: 


^C\G c tanh 



= o^eff 


(B10) 


as given in eq. (50). So it is the specific function f c rc (x, y) 
of eq. (16) that allows us to order the lowest orders in 
the perturbation, on the sole assumption that 0 < $i < 1 , 
and then identifying which value of $i correspond to each 
regime in a e fi. The predictions for the disordered HL 
model are gathered in eq. (50), and we can recover their 
counterparts for the standard HL model by replacing all 
the moments (crj) by crj, as listed in eq. (17). 

Actually, in order to obtain the derivation of the 
Herschel-Bulkley behavior of cr M at a e s < <a c , we need 
to compute the second lowest order of Dhl- We thus 
start from the ansatz: 


Dhlt = x 2 A ) Ci G 0 7t 


1 + G 2 (G 0 7 t ) 1/2 


(BH) 


as suggested by eq. (Bl) for the standard HL model, and 
the same procedure as before leads to the following rela¬ 
tion between C 2 and C\\ 


c 2 = Vcl 


+ Gi (tanh (^) ) - | (a c tanh 2 (^g) ) 
~ (Gicr c tanh (.7)) - | (cr 2 tanh 2 (5^)) 


(B12) 

The resulting predictions for the mean stress gm , and 
specifically for the prefactor A of the stress contribution 
in (G 0 7 t) 1/2 , are given explicitly in Appendix E. 


The expression for C\ — Ci(a e ff) and consequently for 
C 2 = C 2 (Ci) can be considerably simplified in the two 
limiting cases a e n < a c and a e n <a c , and they lead to 
eq. (52). The argument is the following: first, for each 
coupling parameter a e s below a c , we can define a typ¬ 
ical value cr* = 2C±. Then, on the one hand, close to 
a c we have C\ —y oo and cr* —)> oo, so we can safely ne¬ 
glect the contributions of cr c > cr*. For the contributions 
of cr c < cr*, we can approximate the hyperbolic tangent 
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with its Taylor expansion at 2a c /Ci <C 1. On the other 
hand, with a e fi close to zero, we have C\ 0 and cr* —)> 0, 
so we can neglect the contributions of cr c < cr* and use for 
a c > cr* the approximation tanh(2<r c /Ci) ~ 1. In prac¬ 
tice, we can decompose the average (O) in eq. (51) into 
two separate averages, restricted on the contributions 
from a c ^ cr*. So for a e s < a c we have (O) Gc>a * ~ 0 
and (O) ~ (0) a <cr *, whereas at a e ff <C ol c we have 

( °)o c <ai ~ 0 and (°) ~ (°)a c >a* c - These approxima- 
tions eventually lead to the following expressions, on the 
one hand at a e ff < a c : 


C i 

G 2 ~- 

and at a e s <C a c : 


{<*$) 


1/2 


24 (a c - a e fi) 

(^c) 1/4 (°~c) 

2 V 4 X 3 1 / 4 (a c - a e ff) 5/4 


^ <^eff n l <^eff 

(°c) V <*c> 


1/2 


(B13) 


(B14) 


Note finally that the expansions of f Gc (x,Go 7 t/x 2 ) 
given in eqs. (B5)-(B7)-(B6), before averaging over the 
values of cr c , allow us to obtain the predictions for p st (cr c ) 
discussed in sect. IV E. Indeed, we have derived the ex¬ 
pressions listed in eqs. (61)-(68)-(71) by substituting into 
these expansions of f Gc the low-shear-rate diffusion coef¬ 
ficient. 


Appendix C: Normalization condition for a generic 
diffusion coefficient D(cr c ) 


In sect. (IV A), we have derived the normalization con¬ 
dition for the stationary PDF V st (cr c , cr) with the restric¬ 
tion that the diffusion coefficient does not depend on the 
local yield stress cr c , but is rather a global quantity con¬ 
trolling the evolution of the PDF according to eq. (27). 

If the diffusion coefficient is more generically of the 
form D(cr c ,t), in the stationary case the normalization 
condition (40) is: 


Pst (cr c ) 


r s t (&c)r 


fa c 



Gp'yr 

D(ct c )t 


D(a c )r 


(Cl) 


where f Gc is exactly the same function as in eq. (12), for 
instance the parabola (14) in absence of shear rate and 
the function (16) in presence of a constant shear rate. 
Using again the relation (34), we obtain the counterpart 


of eq. (12) for the global plastic activity: 




Gp'yr 

D(<j c )t 


D(a c )r 


= 1 


(C2) 


but this expression does not simplify into eq. (42)-(43), 
and thus the closure relation (44) is modified by the cr c - 
dependence of D(a c ). The previous relation (C2) can be 
used to compute r st , at least numerically if not analyti¬ 
cally, for any choice of p(cr c ) and D(cr c ). 

Nevertheless, for the sake of completeness, we can 
parametrize the stationary diffusion coefficient according 
to: 


D(a c ) = Dd(a c ), / do c d(cr c ) = 1 (C3) 

Jo 

where on the one hand, D is the diffusion coefficient in¬ 
tegrated over all the possible values of cr c (on which we 
could for instance impose a closure relation for r st (D)), 
and on the other hand, d(cr c ) characterizes how the dif¬ 
fusion affects the sites with different values of the local 
yield stress cr c . If such a parametrization is relevant for a 
given amorphous system, then eq. (C2) simply becomes: 


r stT 


f‘OC> 

/ da c p(a c ) 
Jo 


L 777K) 

Dt d(a c ) 


= 1. 


(C4) 

So, combined with the closure relation (36), this last re¬ 
lation provides us with the generalized counterpart of 
eq. (13): 


fa c ( x\jd(a c ),y/d(a c ) 
d(a c ) 




(C5) 


with x = V Dr and y = Go'fr/x 2 . This defines an ‘effec¬ 
tive’ function f e ff y\d(cr c )^ similarly to eq. (44). The 

diffusion coefficient D = Dhl can then be determined 
uniquely as a function of the shear rate 7 and the ef¬ 
fective coupling parameter a e s. Note at last that the 
shape of d(a c ) should be justified separately, as it is here 
introduced as an arbitrary input of the model. 


Appendix D: Explicit analytical expressions for the 
stationary case at fixed diffusion coefficient 

In this section we sketch the derivation and give the 
explicit expressions of the stationary solution of the dis¬ 
ordered HL evolution equation (27), on the one hand the 
complete PDFs and on the other hand the corresponding 
mean stress, at fixed diffusion coefficient (in other words, 
before using any specific closure relation for r Bt m- 
The equation of the stationary joint PDF decomposes 
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into the following structure, respectively on |cr| > cr c and 

M < cr c : 

[dl - p 0 d a \ P(a) = 0 => P{a) = a e /3o<T + c 2 
[dl - p 0 d a \ P(a) = ® =*► P{a) = h e /3(_)<T + c 2 


+ -j^p and the constants 

{ci,ci,C 2 } fixed by the boundary conditions at 
a G {—oo,—cr c ,0,cr c ,oo}. Adapting first these solutions 
to our notations with /3 0 = y = , the joint PDF 

V st (cr c ,a) can be decomposed into: 

Pst(cr c ,a) = p{i k((t c )P< r c (v) (Dl) 


with %) = ^ ± a / 


as announced in sect. IV B, with f st (cr c ) = r st p(cr c ) ac¬ 
cording to eq. (34), and 


PaSe) = < 


^(-) e +/3(-)Cr c 


,y(a-a c ) _)_ £(+) 


fl(-) 


t((±!g-/3(+)<T 0 e !/(o'+o-c) _i_ ^(~) 

» L ^(+>. 


e 


£(+) C 


: a > a c 
: 0 < a < a c 

: — <7 C < a < 0 
: a < <7 c 


( D2 ) 

The definition of the partial plastic activity T st (a c ) in 
eq. (30) allows one to determine the normalization factor: 


k(a c ) = Dr 

= [P(+) e 


/ da p a c (a) 

J \(j\>(Jr 

P{~) e 


|cr|>cr c (D3) 

j pP( — ) ct c — Q ~/3 {+) a c I" 1 


The function f ac (x,y) is then defined with respect to 
the dynamical distribution of local yield stress in the 
stationary case, p(a c ) = p(<j c )^r fa c {x, y), whose defi¬ 
nition (32) implies that 

f ac (x,y) = k(cr c ) / dap ac (a) (D4) 

Jr 

which is thus exactly the same expression (16) as for the 
standard HL model. From the global normalization of 
the PDF, we conclude that the global plastic activity at 

fixed D is given by T st r = Dr/ ^fa c (x,y)^, as stated in 

eqs. (42)-(43). Once f Gc (x,y) is known, the dynamical 
distribution of local yield stress p(cr c ) can be fully deter¬ 
mined according to eq. (45). 

The main novelty in eq. (Dl), compared to previ¬ 
ous references on the standard HL model [21, 24], is 
that the global plastic activity at fixed diffusion coef¬ 
ficient, r st (D), is replaced by its partial counterpart 
r st (cr c ) = r st p(cr c ). Moreover, we have explicitly kept 


the ratio r st /D, with the global plastic activity fixed by 
eqs. (42)-(43)-(12); so the solution (Dl) remains valid for 
any closure relation, and in particular for the HL closure 
relation (35). In the latter case, that we have studied 
throughout this paper, the ratio r st /D can simply be 
replaced by l/a e ff- 


Since all the dependences on the local yield stress a c 
have been made explicit, the stress PDF P st (cr) can be 
computed by integrating V s t (<t c , &) over the possible val¬ 
ues of local yield stress. Nevertheless, for an arbitrary 
a priori distribution p(cr c ), no explicit expression can be 
written down, because of the cr c -dependence of the stress 
division itself (|cr I ^<7c). 


We come at last to the prediction for the mean stress 
<jM(x,y), with x 2 = Dr and y = Gq^t/x 2 fixed, using 

D/r st (D) = ^ fa c (x,y ■)) according to eqs. (42)-(43). We 

distinguish the contributions at fixed local yield stress of 
overstressed and understressed regions: 


< 

\ »£ nd,r> 


-r-r — 7 — 7\(k(o'c) ft i ^ da a Pa (a)) 
(/•tcOc.D) \ v ; • P<T = V V 

(D5) 


We start with the contribution of the overstressed re¬ 
gions: 


k(a c ) / daapa c {a) = fa c (x,y)x 2 y (D 6 ) 
J\cr\>cr c 

which, combined to eq. (45), leads to: 

nOO 

^(over) = / da c p(a c ) Go^T GojT (D7) 

Jo 

We emphasize that this result does not depend on a spe¬ 
cific choice for the closure relation r st (D), it stems solely 
from the specific functional of the stationary joint PDF 
P(cr c , a). We turn now to the contribution of the under¬ 
stressed regions: 


(<r c ) / da aPa c {a) 
J |cr| <cr r 


al . a, !-!- 


= ^ + 


(A + V 1 + A^) tanh (T) 


^ y 2 tanh (^) + 1 + 

{jf-hc {x,y) + x 2S j 1 


4 

x 2 y 2 


2a r . 


y v /rT ^¥ +tanh (T) 

(D 8 ) 
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which leads to 


(under) 

°M 


1 

y 


( a c)/ 2 ~ (L c (x,y)^ + X 2 

(f<Tc( x ,y)) 


+ 


y 2 (fa c (x,y) 



2(7, 


1 + A? + tanh (79' 


(D9) 


where we can recognize a c = (cr;?) /2 in the first term. 
The total mean stress can eventually be computed by 
combining eqs. (D7) and (D9) into 


o M 


(over) (under) 

a M + <J M 


(DIO) 


Moreover, while discussing the assumption of a typical 
fixed rate 1/r in sect. V A, we have suggested the al¬ 
ternative definition of the ‘macroscopic’ stress given in 
eq. (82). It simply consists in the replacement of crj^ ver ^ 
by f^°d<j c g c T st (ct c )t, and hence at fixed (x,y): 


g 


corr 

M 


(<7 C ) a; 2 

(f* c ( x ,y)) 


, (under) 


(Dll) 


The limit of low shear rate 7 —>• 0 of ctm, with the 
HL closure relation (36), is discussed in the main text 
in sect. IV D. Note finally that, before performing any 
Taylor expansion of gm at small 7 , it is crucial not to 
replace (^f (Tc (x 1 y)^ by ce e ff, in order to capture correctly 

the lowest orders in the perturbation; the HL closure 
relation will in fact already be encoded in the diffusion 
coefficient Dhl itself. 


Appendix E: Herschel-Bulkley behavior in the 
disordered HL model at o e ff < ol c 

As discussed in sect. IV D, for a e ff < a c the mean stress 
displays a Herschel-Bulkley behavior a low shear rate: 

(7m (7 ^ 0) <7y + A (Go7t) 1/2 

The macroscopic yield stress Gy is simply obtained us¬ 
ing the lowest-order expansion of the diffusion coefficient 
(using the minimal ansatz of eq. (B3)), the prefactor A 
involves its second-order expansion (using the ansatz of 
eq. (Bll)). In this appendix, we give explicitly the ex¬ 
pressions of these two parameters of the Herschel-Bulkley 


behavior of exponent 1/2 predicted by the disordered HL 
model for a e ff < ol c . 

We first substitute the ansatz for Dhl given in 
eq. (Bll) into the exact expression for the mean stress 
gm of eq. (DIO), and expand the result at small G 07 T. 
We obtain at O ( 7 0 ): 


*y(Cl) = Ci 


{AA1 

Cl (a c tanh (^-) ) 


(El) 


and at O ( 7 1 / 2 ) • 

— - Cl (tanh ( 5 ^-)} + I (( 7 C tanh 2 (^-) 


A = C 1 ’ 2 ^ 


+ 2G2 


(7 C tanh ( 5 ^)) 

A 1 - C1 ( a c tanh (^-) ) - \ (a 2 c tanh 2 (^-) ) 


, tanh ( - 


(*)) 


(E2) 


with C 2 {Ci) given by eq. (B 12 ), and C\(a e ^) by 
eq. (BIO). 

Secondly, we can simplify this expression in the usual 
limiting cases of a e s < a c and a e s <C <a c , as announced in 
eq. (56). On the one hand, we use eq. (B13) at a e ff < a c‘- 


Gy 


(“) 7 c )/ 7 c > («) ( a c - a eff ) 1/2 <( 7^> 1/2 


12C 


V6 ((7 C 2 ) 

l" 3 / 4 /cr 4 \ 3/4 /, 


A ~ 77^3/2 _ ( a c ~ «eff) 3/4 (7) (<7c) 

~ ((7?) ~ 23/2 x 6 3/ 4 (<7 c 2 ) 


(E3) 


and on the other hand, we use eq. (B14) at a e ff a c 
(“) 7c) / 7c) _ (52) a c - tteff 


<7y 


7c) 


1 Ml 

, 2 ((7 C ) 2 


77: 


\ _ ( 7 ) 9 

, 2 ((7 C ) 2 / \7c)/ 


1/2 


(E4) 


The predictions of the standard HL model, presented 
in sect. Ill, respectively eqs. (20)-(22)-(23), are of course 
recovered by removing all the averages over p(g c ) from 
the three last equations. We emphasize that the absence 
of averages over p(g c ) allows one to simplify considerably 
these expressions, removing in particular the non-trivial 
combinations of moments ( 0 /?). 
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